1 SUBROUTINE DGELSY_F95( A, B, RANK, JPVT, RCOND, INFO ) 2! 3! -- LAPACK95 interface driver routine (version 3.0) -- 4! UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK 5! September, 2000 6! 7! .. USE STATEMENTS .. 8 USE LA_PRECISION, ONLY: WP => DP 9 USE LA_AUXMOD, ONLY: ERINFO 10 USE F77_LAPACK, ONLY: GELSY_F77 => LA_GELSY 11! .. IMPLICIT STATEMENT .. 12 IMPLICIT NONE 13! .. SCALAR ARGUMENTS .. 14 INTEGER, INTENT(OUT), OPTIONAL :: RANK 15 INTEGER, INTENT(OUT), OPTIONAL :: INFO 16 REAL(WP), INTENT(IN), OPTIONAL :: RCOND 17! .. ARRAY ARGUMENTS .. 18 INTEGER, INTENT(INOUT), OPTIONAL, TARGET :: JPVT(:) 19 REAL(WP), INTENT(INOUT) :: A(:,:), B(:,:) 20!---------------------------------------------------------------------- 21! 22! Purpose 23! ======= 24! 25! LA_GELSY computes the minimum-norm least squares solution to one 26! or more real or complex linear systems A*x = b using a complete 27! orthogonal factorization of A. Matrix A is rectangular and may be 28! rankdeficient. The vectors b and corresponding solution vectors x are 29! the columns of matrices denoted B and X, respectively. 30! The routine computes a QR factorization of A with column pivoting: 31! A * P = Q * [ R11 R12 ] 32! [ 0 R22 ] 33! where R11 is the largest leading submatrix whose estimated condition 34! number is less than 1/RCOND. The order of R11, RANK, is the effective 35! rank of A. R22 is considered to be negligible, and R12 is annihilated 36! by orthogonal (unitary) transformations from the right, yielding the 37! complete orthogonal (unitary) factorization 38! A * P = Q * [ T11 0 ] * Z 39! [ 0 0 ] 40! The minimum-norm least squares solution is then 41! x = P * Z^H [ T11^-1 * Q1^H * b ] 42! [ 0 ] 43! where Q1 consists of the first RANK columns of Q. 44! 45! ========= 46! 47! SUBROUTINE LA_GELSY( A, B, RANK=rank, & 48! JPVT= jpvt, RCOND= rcond, INFO= info ) 49! <type>(<wp>), INTENT(INOUT) :: A(:,:), <rhs> 50! INTEGER, INTENT(OUT), OPTIONAL :: RANK 51! INTEGER, INTENT(INOUT), OPTIONAL :: JPVT(:) 52! REAL(<wp>), INTENT(IN), OPTIONAL :: RCOND 53! INTEGER, INTENT(OUT), OPTIONAL :: INFO 54! where 55! <type> ::= REAL | COMPLEX 56! <wp> ::= KIND(1.0) | KIND(1.0D0) 57! <rhs> ::= B(:,:) | B(:) 58! 59! Arguments 60! ========= 61! 62! A (input/output) REAL or COMPLEX array, shape (:,:). 63! On entry, the matrix A. 64! On exit, A has been overwritten by details of its complete 65! orthogonal factorization. 66! B (input/output) REAL or COMPLEX array, shape (:,:) with 67! size(B,1) = max(size(A,1),size(A,2)) or shape (:) with 68! size(B) = max(size(A,1), size(A,2)). 69! On entry, the matrix B. 70! On exit, rows 1 to size(A,2) contain the solution matrix X . 71! If size(A,1) >= size(A,2) and RANK = size(A,2), the residual 72! sum-of-squares for the solution vector in a column of B is 73! given by the sum of squares of elements in rows size(A,2)+1 : 74! size(A,1) of that column. 75! RANK Optional (output) INTEGER. 76! The effective rank of A, i.e., the order of the submatrix R11. 77! This is the same as the order of the submatrix T11 in the 78! complete orthogonal factorization of A. 79! JPVT Optional (input/output) INTEGER array, shape (:) with 80! size(JPVT) = size(A,2). 81! On entry, if JPVT(i) /= 0, the i-th column of A is an initial 82! column, otherwise it is a free column. 83! Before the QR factorization of A, all initial columns are 84! permuted to the leading positions; only the remaining free 85! columns are moved as a result of column pivoting during the 86! factorization. 87! On exit, if JPVT(i) = k, then the i-th column of the matrix 88! product A*P was the k-th column of A. 89! RCOND Optional (input) REAL. 90! RCOND is used to determine the effective rank of A. This is 91! defined as the order of the largest leading triangular 92! submatrix R11 in the QR factorization of A, with pivoting, 93! whose estimated condition number < 1/RCOND. 94! Default value: 10*max(size(A,1),size(A,2))*BEPSILON(1.0_<wp>), 95! where <wp> is the working precision. 96! INFO Optional (output) INTEGER. 97! = 0: successful exit 98! < 0: if INFO = -i, the i-th argument had an illegal value 99! If INFO is not present and an error occurs, then the program 100! is terminated with an error message. 101!---------------------------------------------------------------------- 102! .. PARAMETERS .. 103 CHARACTER(LEN=8), PARAMETER :: SRNAME = 'LA_GELSY' 104! .. LOCAL SCALARS .. 105 INTEGER :: LINFO, ISTAT, ISTAT1, LWORK, N, M, MN, NRHS, LRANK, SJPVT 106 REAL(WP) :: LRCOND 107! .. LOCAL POINTERS .. 108 INTEGER, POINTER :: LJPVT(:) 109 REAL(WP), POINTER :: WORK(:) 110 REAL(WP) :: WORKMIN(1) 111! .. INTRINSIC FUNCTIONS .. 112 INTRINSIC SIZE, PRESENT, MAX, MIN, EPSILON 113! .. EXECUTABLE STATEMENTS .. 114 LINFO = 0; ISTAT = 0; M = SIZE(A,1); N = SIZE(A,2); NRHS = SIZE(B,2) 115 MN = MIN(M,N) 116 IF( PRESENT(RCOND) )THEN; LRCOND = RCOND; ELSE 117 LRCOND = 100*EPSILON(1.0_WP) ; ENDIF 118 IF( PRESENT(JPVT) )THEN; SJPVT = SIZE(JPVT); ELSE; SJPVT = N; ENDIF 119! .. TEST THE ARGUMENTS 120 IF( M < 0 .OR. N < 0 ) THEN; LINFO = -1 121 ELSE IF( SIZE( B, 1 ) /= MAX(1,M,N) .OR. NRHS < 0 ) THEN; LINFO = -2 122 ELSE IF( SJPVT /= N ) THEN; LINFO = -4 123 ELSE IF( LRCOND <= 0.0_WP ) THEN; LINFO = -5 124 ELSE 125 IF( PRESENT(JPVT) )THEN; LJPVT => JPVT 126 ELSE; ALLOCATE( LJPVT(N), STAT = ISTAT ); LJPVT = 0; END IF 127 128 129! .. DETERMINE THE WORKSPACE .. 130! .. QUERING THE SIZE OF WORKSPACE .. 131 LWORK = -1 132 CALL GELSY_F77( M, N, NRHS, A, MAX(1,M), B, MAX(1,M,N), & 133 & LJPVT, LRCOND, LRANK, WORKMIN, LWORK, LINFO ) 134 LWORK = WORKMIN(1) 135 IF( ISTAT == 0 ) THEN 136 ALLOCATE( WORK(LWORK), STAT = ISTAT ) 137 IF( ISTAT /= 0 ) CALL ERINFO( -200, SRNAME, LINFO ) 138 END IF 139 140 IF ( ISTAT == 0 ) THEN 141 CALL GELSY_F77( M, N, NRHS, A, MAX(1,M), B, MAX(1,M,N), & 142 & LJPVT, LRCOND, LRANK, WORK, LWORK, LINFO ) 143 ELSE; LINFO = -100; END IF 144 IF( PRESENT(RANK) ) RANK = LRANK 145 IF( PRESENT(JPVT) ) JPVT = LJPVT 146 DEALLOCATE(WORK, STAT = ISTAT1 ) 147 END IF 148 CALL ERINFO( LINFO, SRNAME, INFO, ISTAT ) 149 END SUBROUTINE DGELSY_F95 150