Purpose ======= LA_GELSY computes the minimum-norm least squares solution to one or more real or complex linear systems A*x = b using a complete orthogonal factorization of A. Matrix A is rectangular and may be rankdeficient. The vectors b and corresponding solution vectors x are the columns of matrices denoted B and X, respectively. The routine computes a QR factorization of A with column pivoting: A * P = Q * [ R11 R12 ] [ 0 R22 ] where R11 is the largest leading submatrix whose estimated condition number is less than 1/RCOND. The order of R11, RANK, is the effective rank of A. R22 is considered to be negligible, and R12 is annihilated by orthogonal (unitary) transformations from the right, yielding the complete orthogonal (unitary) factorization A * P = Q * [ T11 0 ] * Z [ 0 0 ] The minimum-norm least squares solution is then x = P * Z^H [ T11^-1 * Q1^H * b ] [ 0 ] where Q1 consists of the first RANK columns of Q. ========= SUBROUTINE LA_GELSY( A, B, RANK=rank, & JPVT= jpvt, RCOND= rcond, INFO= info ) (), INTENT(INOUT) :: A(:,:), INTEGER, INTENT(OUT), OPTIONAL :: RANK INTEGER, INTENT(INOUT), OPTIONAL :: JPVT(:) REAL(), INTENT(IN), OPTIONAL :: RCOND INTEGER, INTENT(OUT), OPTIONAL :: INFO where ::= REAL | COMPLEX ::= KIND(1.0) | KIND(1.0D0) ::= B(:,:) | B(:) Arguments ========= A (input/output) REAL or COMPLEX array, shape (:,:). On entry, the matrix A. On exit, A has been overwritten by details of its complete orthogonal factorization. B (input/output) REAL or COMPLEX array, shape (:,:) with size(B,1) = max(size(A,1),size(A,2)) or shape (:) with size(B) = max(size(A,1), size(A,2)). On entry, the matrix B. On exit, rows 1 to size(A,2) contain the solution matrix X . If size(A,1) >= size(A,2) and RANK = size(A,2), the residual sum-of-squares for the solution vector in a column of B is given by the sum of squares of elements in rows size(A,2)+1 : size(A,1) of that column. RANK Optional (output) INTEGER. The effective rank of A, i.e., the order of the submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of A. JPVT Optional (input/output) INTEGER array, shape (:) with size(JPVT) = size(A,2). On entry, if JPVT(i) /= 0, the i-th column of A is an initial column, otherwise it is a free column. Before the QR factorization of A, all initial columns are permuted to the leading positions; only the remaining free columns are moved as a result of column pivoting during the factorization. On exit, if JPVT(i) = k, then the i-th column of the matrix product A*P was the k-th column of A. RCOND Optional (input) REAL. RCOND is used to determine the effective rank of A. This is defined as the order of the largest leading triangular submatrix R11 in the QR factorization of A, with pivoting, whose estimated condition number < 1/RCOND. Default value: 10*max(size(A,1),size(A,2))*BEPSILON(1.0_), where is the working precision. INFO Optional (output) INTEGER. = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value If INFO is not present and an error occurs, then the program is terminated with an error message.