DULSIA(3)

SYNOPSIS
SUBROUTINE DULSIA (A, MDA, M, N, B, MDB, NB, RE, AE, KEY, MODE, + NP, KRANK, KSURE, RNORM, W, LW, IWORK, LIW, INFO)

NAME
dulsia - Solve an underdetermined linear system of equations by performing an LQ factorization of the matrix using Householder transformations. Emphasis is put on detecting possible rank deficiency.

LIBRARY
SLATEC

CATEGORY
D9

TYPE
DOUBLE PRECISION (ULSIA-S, DULSIA-D)

KEYWORDS
LINEAR LEAST SQUARES, LQ FACTORIZATION, UNDERDETERMINED LINEAR SYSTEM

AUTHOR
Manteuffel, T. A., (LANL)

DESCRIPTION
DULSIA computes the minimal length solution(s) to the problem AX=B where A is an M by N matrix with M.LE.N and B is the M by NB matrix of right hand sides. User input bounds on the uncertainty in the elements of A are used to detect numerical rank deficiency. The algorithm employs a row and column pivot strategy to minimize the growth of uncertainty and round-off errors. DULSIA requires (MDA+1)*N + (MDB+1)*NB + 6*M dimensioned space ****************************************************************** * * * WARNING - All input arrays are changed on exit. * * * ****************************************************************** Input.. All TYPE REAL variables are DOUBLE PRECISION A(,) Linear coefficient matrix of AX=B, with MDA the MDA,M,N actual first dimension of A in the calling program. M is the row dimension (no. of EQUATIONS of the problem) and N the col dimension (no. of UNKNOWNS). Must have MDA.GE.M and M.LE.N. B(,) Right hand side(s), with MDB the actual first MDB,NB dimension of B in the calling program. NB is the number of M by 1 right hand sides. Since the solution is returned in B, must have MDB.GE.N. If NB = 0, B is never accessed. ****************************************************************** * * * Note - Use of RE and AE are what make this * * code significantly different from * * other linear least squares solvers. * * However, the inexperienced user is * * advised to set RE=0.,AE=0.,KEY=0. * * * ****************************************************************** RE(),AE(),KEY RE() RE() is a vector of length N such that RE(I) is the maximum relative uncertainty in row I of the matrix A. The values of RE() must be between 0 and 1. A minimum of 10*machine precision will be enforced. AE() AE() is a vector of length N such that AE(I) is the maximum absolute uncertainty in row I of the matrix A. The values of AE() must be greater than or equal to 0. KEY For ease of use, RE and AE may be input as either vectors or scalars. If a scalar is input, the algo- rithm will use that value for each column of A. The parameter KEY indicates whether scalars or vectors are being input. KEY=0 RE scalar AE scalar KEY=1 RE vector AE scalar KEY=2 RE scalar AE vector KEY=3 RE vector AE vector MODE The integer MODE indicates how the routine is to react if rank deficiency is detected. If MODE = 0 return immediately, no solution 1 compute truncated solution 2 compute minimal length least squares sol The inexperienced user is advised to set MODE=0 NP The first NP rows of A will not be interchanged with other rows even though the pivot strategy would suggest otherwise. The inexperienced user is advised to set NP=0. WORK() A real work array dimensioned 5*M. However, if RE or AE have been specified as vectors, dimension WORK 4*M. If both RE and AE have been specified as vectors, dimension WORK 3*M. LW Actual dimension of WORK IWORK() Integer work array dimensioned at least N+M. LIW Actual dimension of IWORK. INFO Is a flag which provides for the efficient solution of subsequent problems involving the same A but different B. If INFO = 0 original call INFO = 1 subsequent calls On subsequent calls, the user must supply A, KRANK, LW, IWORK, LIW, and the first 2*M locations of WORK as output by the original call to DULSIA. MODE must be equal to the value of MODE in the original call. If MODE.LT.2, only the first N locations of WORK are accessed. AE, RE, KEY, and NP are not accessed. Output..All TYPE REAL variables are DOUBLE PRECISION A(,) Contains the lower triangular part of the reduced matrix and the transformation information. It togeth with the first M elements of WORK (see below) completely specify the LQ factorization of A. B(,) Contains the N by NB solution matrix for X. KRANK,KSURE The numerical rank of A, based upon the relative and absolute bounds on uncertainty, is bounded above by KRANK and below by KSURE. The algorithm returns a solution based on KRANK. KSURE provides an indication of the precision of the rank. RNORM() Contains the Euclidean length of the NB residual vectors B(I)-AX(I), I=1,NB. If the matrix A is of full rank, then RNORM=0.0. WORK() The first M locations of WORK contain values necessary to reproduce the Householder transformation. IWORK() The first N locations contain the order in which the columns of A were used. The next M locations contain the order in which the rows of A were used. INFO Flag to indicate status of computation on completion -1 Parameter error(s) 0 - Rank deficient, no solution 1 - Rank deficient, truncated solution 2 - Rank deficient, minimal length least squares sol 3 - Numerical rank 0, zero solution 4 - Rank .LT. NP 5 - Full rank

REFERENCES
T. Manteuffel, An interval analysis approach to rank determination in linear least squares problems, Report SAND80-0655, Sandia Laboratories, June 1980.

ROUTINES CALLED
D1MACH, DU11US, DU12US, XERMSG

REVISION HISTORY (YYMMDD)
810801 DATE WRITTEN
890831 Modified array declarations. (WRB)
891006 Cosmetic changes to prologue. (WRB)
891009 Removed unreferenced variable. (WRB)
891009 REVISION DATE from Version 3.2
891214 Prologue converted to Version 4.0 format. (BAB)
900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
900510 Fixed an error message. (RWC)
920501 Reformatted the REFERENCES section. (WRB)