root / src / lapack / double / dgelsd.f @ 1
Historique | Voir | Annoter | Télécharger (18,5 ko)
| 1 | 1 | equemene | SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, |
|---|---|---|---|
| 2 | 1 | equemene | $ WORK, LWORK, IWORK, INFO ) |
| 3 | 1 | equemene | * |
| 4 | 1 | equemene | * -- LAPACK driver routine (version 3.2.2) -- |
| 5 | 1 | equemene | * -- LAPACK is a software package provided by Univ. of Tennessee, -- |
| 6 | 1 | equemene | * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
| 7 | 1 | equemene | * June 2010 |
| 8 | 1 | equemene | * |
| 9 | 1 | equemene | * .. Scalar Arguments .. |
| 10 | 1 | equemene | INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK |
| 11 | 1 | equemene | DOUBLE PRECISION RCOND |
| 12 | 1 | equemene | * .. |
| 13 | 1 | equemene | * .. Array Arguments .. |
| 14 | 1 | equemene | INTEGER IWORK( * ) |
| 15 | 1 | equemene | DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) |
| 16 | 1 | equemene | * .. |
| 17 | 1 | equemene | * |
| 18 | 1 | equemene | * Purpose |
| 19 | 1 | equemene | * ======= |
| 20 | 1 | equemene | * |
| 21 | 1 | equemene | * DGELSD computes the minimum-norm solution to a real linear least |
| 22 | 1 | equemene | * squares problem: |
| 23 | 1 | equemene | * minimize 2-norm(| b - A*x |) |
| 24 | 1 | equemene | * using the singular value decomposition (SVD) of A. A is an M-by-N |
| 25 | 1 | equemene | * matrix which may be rank-deficient. |
| 26 | 1 | equemene | * |
| 27 | 1 | equemene | * Several right hand side vectors b and solution vectors x can be |
| 28 | 1 | equemene | * handled in a single call; they are stored as the columns of the |
| 29 | 1 | equemene | * M-by-NRHS right hand side matrix B and the N-by-NRHS solution |
| 30 | 1 | equemene | * matrix X. |
| 31 | 1 | equemene | * |
| 32 | 1 | equemene | * The problem is solved in three steps: |
| 33 | 1 | equemene | * (1) Reduce the coefficient matrix A to bidiagonal form with |
| 34 | 1 | equemene | * Householder transformations, reducing the original problem |
| 35 | 1 | equemene | * into a "bidiagonal least squares problem" (BLS) |
| 36 | 1 | equemene | * (2) Solve the BLS using a divide and conquer approach. |
| 37 | 1 | equemene | * (3) Apply back all the Householder tranformations to solve |
| 38 | 1 | equemene | * the original least squares problem. |
| 39 | 1 | equemene | * |
| 40 | 1 | equemene | * The effective rank of A is determined by treating as zero those |
| 41 | 1 | equemene | * singular values which are less than RCOND times the largest singular |
| 42 | 1 | equemene | * value. |
| 43 | 1 | equemene | * |
| 44 | 1 | equemene | * The divide and conquer algorithm makes very mild assumptions about |
| 45 | 1 | equemene | * floating point arithmetic. It will work on machines with a guard |
| 46 | 1 | equemene | * digit in add/subtract, or on those binary machines without guard |
| 47 | 1 | equemene | * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or |
| 48 | 1 | equemene | * Cray-2. It could conceivably fail on hexadecimal or decimal machines |
| 49 | 1 | equemene | * without guard digits, but we know of none. |
| 50 | 1 | equemene | * |
| 51 | 1 | equemene | * Arguments |
| 52 | 1 | equemene | * ========= |
| 53 | 1 | equemene | * |
| 54 | 1 | equemene | * M (input) INTEGER |
| 55 | 1 | equemene | * The number of rows of A. M >= 0. |
| 56 | 1 | equemene | * |
| 57 | 1 | equemene | * N (input) INTEGER |
| 58 | 1 | equemene | * The number of columns of A. N >= 0. |
| 59 | 1 | equemene | * |
| 60 | 1 | equemene | * NRHS (input) INTEGER |
| 61 | 1 | equemene | * The number of right hand sides, i.e., the number of columns |
| 62 | 1 | equemene | * of the matrices B and X. NRHS >= 0. |
| 63 | 1 | equemene | * |
| 64 | 1 | equemene | * A (input) DOUBLE PRECISION array, dimension (LDA,N) |
| 65 | 1 | equemene | * On entry, the M-by-N matrix A. |
| 66 | 1 | equemene | * On exit, A has been destroyed. |
| 67 | 1 | equemene | * |
| 68 | 1 | equemene | * LDA (input) INTEGER |
| 69 | 1 | equemene | * The leading dimension of the array A. LDA >= max(1,M). |
| 70 | 1 | equemene | * |
| 71 | 1 | equemene | * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) |
| 72 | 1 | equemene | * On entry, the M-by-NRHS right hand side matrix B. |
| 73 | 1 | equemene | * On exit, B is overwritten by the N-by-NRHS solution |
| 74 | 1 | equemene | * matrix X. If m >= n and RANK = n, the residual |
| 75 | 1 | equemene | * sum-of-squares for the solution in the i-th column is given |
| 76 | 1 | equemene | * by the sum of squares of elements n+1:m in that column. |
| 77 | 1 | equemene | * |
| 78 | 1 | equemene | * LDB (input) INTEGER |
| 79 | 1 | equemene | * The leading dimension of the array B. LDB >= max(1,max(M,N)). |
| 80 | 1 | equemene | * |
| 81 | 1 | equemene | * S (output) DOUBLE PRECISION array, dimension (min(M,N)) |
| 82 | 1 | equemene | * The singular values of A in decreasing order. |
| 83 | 1 | equemene | * The condition number of A in the 2-norm = S(1)/S(min(m,n)). |
| 84 | 1 | equemene | * |
| 85 | 1 | equemene | * RCOND (input) DOUBLE PRECISION |
| 86 | 1 | equemene | * RCOND is used to determine the effective rank of A. |
| 87 | 1 | equemene | * Singular values S(i) <= RCOND*S(1) are treated as zero. |
| 88 | 1 | equemene | * If RCOND < 0, machine precision is used instead. |
| 89 | 1 | equemene | * |
| 90 | 1 | equemene | * RANK (output) INTEGER |
| 91 | 1 | equemene | * The effective rank of A, i.e., the number of singular values |
| 92 | 1 | equemene | * which are greater than RCOND*S(1). |
| 93 | 1 | equemene | * |
| 94 | 1 | equemene | * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) |
| 95 | 1 | equemene | * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
| 96 | 1 | equemene | * |
| 97 | 1 | equemene | * LWORK (input) INTEGER |
| 98 | 1 | equemene | * The dimension of the array WORK. LWORK must be at least 1. |
| 99 | 1 | equemene | * The exact minimum amount of workspace needed depends on M, |
| 100 | 1 | equemene | * N and NRHS. As long as LWORK is at least |
| 101 | 1 | equemene | * 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2, |
| 102 | 1 | equemene | * if M is greater than or equal to N or |
| 103 | 1 | equemene | * 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2, |
| 104 | 1 | equemene | * if M is less than N, the code will execute correctly. |
| 105 | 1 | equemene | * SMLSIZ is returned by ILAENV and is equal to the maximum |
| 106 | 1 | equemene | * size of the subproblems at the bottom of the computation |
| 107 | 1 | equemene | * tree (usually about 25), and |
| 108 | 1 | equemene | * NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) |
| 109 | 1 | equemene | * For good performance, LWORK should generally be larger. |
| 110 | 1 | equemene | * |
| 111 | 1 | equemene | * If LWORK = -1, then a workspace query is assumed; the routine |
| 112 | 1 | equemene | * only calculates the optimal size of the WORK array, returns |
| 113 | 1 | equemene | * this value as the first entry of the WORK array, and no error |
| 114 | 1 | equemene | * message related to LWORK is issued by XERBLA. |
| 115 | 1 | equemene | * |
| 116 | 1 | equemene | * IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK)) |
| 117 | 1 | equemene | * LIWORK >= max(1, 3 * MINMN * NLVL + 11 * MINMN), |
| 118 | 1 | equemene | * where MINMN = MIN( M,N ). |
| 119 | 1 | equemene | * On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK. |
| 120 | 1 | equemene | * |
| 121 | 1 | equemene | * INFO (output) INTEGER |
| 122 | 1 | equemene | * = 0: successful exit |
| 123 | 1 | equemene | * < 0: if INFO = -i, the i-th argument had an illegal value. |
| 124 | 1 | equemene | * > 0: the algorithm for computing the SVD failed to converge; |
| 125 | 1 | equemene | * if INFO = i, i off-diagonal elements of an intermediate |
| 126 | 1 | equemene | * bidiagonal form did not converge to zero. |
| 127 | 1 | equemene | * |
| 128 | 1 | equemene | * Further Details |
| 129 | 1 | equemene | * =============== |
| 130 | 1 | equemene | * |
| 131 | 1 | equemene | * Based on contributions by |
| 132 | 1 | equemene | * Ming Gu and Ren-Cang Li, Computer Science Division, University of |
| 133 | 1 | equemene | * California at Berkeley, USA |
| 134 | 1 | equemene | * Osni Marques, LBNL/NERSC, USA |
| 135 | 1 | equemene | * |
| 136 | 1 | equemene | * ===================================================================== |
| 137 | 1 | equemene | * |
| 138 | 1 | equemene | * .. Parameters .. |
| 139 | 1 | equemene | DOUBLE PRECISION ZERO, ONE, TWO |
| 140 | 1 | equemene | PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) |
| 141 | 1 | equemene | * .. |
| 142 | 1 | equemene | * .. Local Scalars .. |
| 143 | 1 | equemene | LOGICAL LQUERY |
| 144 | 1 | equemene | INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ, |
| 145 | 1 | equemene | $ LDWORK, LIWORK, MAXMN, MAXWRK, MINMN, MINWRK, |
| 146 | 1 | equemene | $ MM, MNTHR, NLVL, NWORK, SMLSIZ, WLALSD |
| 147 | 1 | equemene | DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM |
| 148 | 1 | equemene | * .. |
| 149 | 1 | equemene | * .. External Subroutines .. |
| 150 | 1 | equemene | EXTERNAL DGEBRD, DGELQF, DGEQRF, DLABAD, DLACPY, DLALSD, |
| 151 | 1 | equemene | $ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR, XERBLA |
| 152 | 1 | equemene | * .. |
| 153 | 1 | equemene | * .. External Functions .. |
| 154 | 1 | equemene | INTEGER ILAENV |
| 155 | 1 | equemene | DOUBLE PRECISION DLAMCH, DLANGE |
| 156 | 1 | equemene | EXTERNAL ILAENV, DLAMCH, DLANGE |
| 157 | 1 | equemene | * .. |
| 158 | 1 | equemene | * .. Intrinsic Functions .. |
| 159 | 1 | equemene | INTRINSIC DBLE, INT, LOG, MAX, MIN |
| 160 | 1 | equemene | * .. |
| 161 | 1 | equemene | * .. Executable Statements .. |
| 162 | 1 | equemene | * |
| 163 | 1 | equemene | * Test the input arguments. |
| 164 | 1 | equemene | * |
| 165 | 1 | equemene | INFO = 0 |
| 166 | 1 | equemene | MINMN = MIN( M, N ) |
| 167 | 1 | equemene | MAXMN = MAX( M, N ) |
| 168 | 1 | equemene | MNTHR = ILAENV( 6, 'DGELSD', ' ', M, N, NRHS, -1 ) |
| 169 | 1 | equemene | LQUERY = ( LWORK.EQ.-1 ) |
| 170 | 1 | equemene | IF( M.LT.0 ) THEN |
| 171 | 1 | equemene | INFO = -1 |
| 172 | 1 | equemene | ELSE IF( N.LT.0 ) THEN |
| 173 | 1 | equemene | INFO = -2 |
| 174 | 1 | equemene | ELSE IF( NRHS.LT.0 ) THEN |
| 175 | 1 | equemene | INFO = -3 |
| 176 | 1 | equemene | ELSE IF( LDA.LT.MAX( 1, M ) ) THEN |
| 177 | 1 | equemene | INFO = -5 |
| 178 | 1 | equemene | ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN |
| 179 | 1 | equemene | INFO = -7 |
| 180 | 1 | equemene | END IF |
| 181 | 1 | equemene | * |
| 182 | 1 | equemene | SMLSIZ = ILAENV( 9, 'DGELSD', ' ', 0, 0, 0, 0 ) |
| 183 | 1 | equemene | * |
| 184 | 1 | equemene | * Compute workspace. |
| 185 | 1 | equemene | * (Note: Comments in the code beginning "Workspace:" describe the |
| 186 | 1 | equemene | * minimal amount of workspace needed at that point in the code, |
| 187 | 1 | equemene | * as well as the preferred amount for good performance. |
| 188 | 1 | equemene | * NB refers to the optimal block size for the immediately |
| 189 | 1 | equemene | * following subroutine, as returned by ILAENV.) |
| 190 | 1 | equemene | * |
| 191 | 1 | equemene | MINWRK = 1 |
| 192 | 1 | equemene | LIWORK = 1 |
| 193 | 1 | equemene | MINMN = MAX( 1, MINMN ) |
| 194 | 1 | equemene | NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ+1 ) ) / |
| 195 | 1 | equemene | $ LOG( TWO ) ) + 1, 0 ) |
| 196 | 1 | equemene | * |
| 197 | 1 | equemene | IF( INFO.EQ.0 ) THEN |
| 198 | 1 | equemene | MAXWRK = 0 |
| 199 | 1 | equemene | LIWORK = 3*MINMN*NLVL + 11*MINMN |
| 200 | 1 | equemene | MM = M |
| 201 | 1 | equemene | IF( M.GE.N .AND. M.GE.MNTHR ) THEN |
| 202 | 1 | equemene | * |
| 203 | 1 | equemene | * Path 1a - overdetermined, with many more rows than columns. |
| 204 | 1 | equemene | * |
| 205 | 1 | equemene | MM = N |
| 206 | 1 | equemene | MAXWRK = MAX( MAXWRK, N+N*ILAENV( 1, 'DGEQRF', ' ', M, N, |
| 207 | 1 | equemene | $ -1, -1 ) ) |
| 208 | 1 | equemene | MAXWRK = MAX( MAXWRK, N+NRHS* |
| 209 | 1 | equemene | $ ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N, -1 ) ) |
| 210 | 1 | equemene | END IF |
| 211 | 1 | equemene | IF( M.GE.N ) THEN |
| 212 | 1 | equemene | * |
| 213 | 1 | equemene | * Path 1 - overdetermined or exactly determined. |
| 214 | 1 | equemene | * |
| 215 | 1 | equemene | MAXWRK = MAX( MAXWRK, 3*N+( MM+N )* |
| 216 | 1 | equemene | $ ILAENV( 1, 'DGEBRD', ' ', MM, N, -1, -1 ) ) |
| 217 | 1 | equemene | MAXWRK = MAX( MAXWRK, 3*N+NRHS* |
| 218 | 1 | equemene | $ ILAENV( 1, 'DORMBR', 'QLT', MM, NRHS, N, -1 ) ) |
| 219 | 1 | equemene | MAXWRK = MAX( MAXWRK, 3*N+( N-1 )* |
| 220 | 1 | equemene | $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, N, -1 ) ) |
| 221 | 1 | equemene | WLALSD = 9*N+2*N*SMLSIZ+8*N*NLVL+N*NRHS+(SMLSIZ+1)**2 |
| 222 | 1 | equemene | MAXWRK = MAX( MAXWRK, 3*N+WLALSD ) |
| 223 | 1 | equemene | MINWRK = MAX( 3*N+MM, 3*N+NRHS, 3*N+WLALSD ) |
| 224 | 1 | equemene | END IF |
| 225 | 1 | equemene | IF( N.GT.M ) THEN |
| 226 | 1 | equemene | WLALSD = 9*M+2*M*SMLSIZ+8*M*NLVL+M*NRHS+(SMLSIZ+1)**2 |
| 227 | 1 | equemene | IF( N.GE.MNTHR ) THEN |
| 228 | 1 | equemene | * |
| 229 | 1 | equemene | * Path 2a - underdetermined, with many more columns |
| 230 | 1 | equemene | * than rows. |
| 231 | 1 | equemene | * |
| 232 | 1 | equemene | MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) |
| 233 | 1 | equemene | MAXWRK = MAX( MAXWRK, M*M+4*M+2*M* |
| 234 | 1 | equemene | $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) |
| 235 | 1 | equemene | MAXWRK = MAX( MAXWRK, M*M+4*M+NRHS* |
| 236 | 1 | equemene | $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, M, -1 ) ) |
| 237 | 1 | equemene | MAXWRK = MAX( MAXWRK, M*M+4*M+( M-1 )* |
| 238 | 1 | equemene | $ ILAENV( 1, 'DORMBR', 'PLN', M, NRHS, M, -1 ) ) |
| 239 | 1 | equemene | IF( NRHS.GT.1 ) THEN |
| 240 | 1 | equemene | MAXWRK = MAX( MAXWRK, M*M+M+M*NRHS ) |
| 241 | 1 | equemene | ELSE |
| 242 | 1 | equemene | MAXWRK = MAX( MAXWRK, M*M+2*M ) |
| 243 | 1 | equemene | END IF |
| 244 | 1 | equemene | MAXWRK = MAX( MAXWRK, M+NRHS* |
| 245 | 1 | equemene | $ ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M, -1 ) ) |
| 246 | 1 | equemene | MAXWRK = MAX( MAXWRK, M*M+4*M+WLALSD ) |
| 247 | 1 | equemene | ! XXX: Ensure the Path 2a case below is triggered. The workspace |
| 248 | 1 | equemene | ! calculation should use queries for all routines eventually. |
| 249 | 1 | equemene | MAXWRK = MAX( MAXWRK, |
| 250 | 1 | equemene | $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) ) |
| 251 | 1 | equemene | ELSE |
| 252 | 1 | equemene | * |
| 253 | 1 | equemene | * Path 2 - remaining underdetermined cases. |
| 254 | 1 | equemene | * |
| 255 | 1 | equemene | MAXWRK = 3*M + ( N+M )*ILAENV( 1, 'DGEBRD', ' ', M, N, |
| 256 | 1 | equemene | $ -1, -1 ) |
| 257 | 1 | equemene | MAXWRK = MAX( MAXWRK, 3*M+NRHS* |
| 258 | 1 | equemene | $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, N, -1 ) ) |
| 259 | 1 | equemene | MAXWRK = MAX( MAXWRK, 3*M+M* |
| 260 | 1 | equemene | $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, M, -1 ) ) |
| 261 | 1 | equemene | MAXWRK = MAX( MAXWRK, 3*M+WLALSD ) |
| 262 | 1 | equemene | END IF |
| 263 | 1 | equemene | MINWRK = MAX( 3*M+NRHS, 3*M+M, 3*M+WLALSD ) |
| 264 | 1 | equemene | END IF |
| 265 | 1 | equemene | MINWRK = MIN( MINWRK, MAXWRK ) |
| 266 | 1 | equemene | WORK( 1 ) = MAXWRK |
| 267 | 1 | equemene | IWORK( 1 ) = LIWORK |
| 268 | 1 | equemene | |
| 269 | 1 | equemene | IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN |
| 270 | 1 | equemene | INFO = -12 |
| 271 | 1 | equemene | END IF |
| 272 | 1 | equemene | END IF |
| 273 | 1 | equemene | * |
| 274 | 1 | equemene | IF( INFO.NE.0 ) THEN |
| 275 | 1 | equemene | CALL XERBLA( 'DGELSD', -INFO ) |
| 276 | 1 | equemene | RETURN |
| 277 | 1 | equemene | ELSE IF( LQUERY ) THEN |
| 278 | 1 | equemene | GO TO 10 |
| 279 | 1 | equemene | END IF |
| 280 | 1 | equemene | * |
| 281 | 1 | equemene | * Quick return if possible. |
| 282 | 1 | equemene | * |
| 283 | 1 | equemene | IF( M.EQ.0 .OR. N.EQ.0 ) THEN |
| 284 | 1 | equemene | RANK = 0 |
| 285 | 1 | equemene | RETURN |
| 286 | 1 | equemene | END IF |
| 287 | 1 | equemene | * |
| 288 | 1 | equemene | * Get machine parameters. |
| 289 | 1 | equemene | * |
| 290 | 1 | equemene | EPS = DLAMCH( 'P' ) |
| 291 | 1 | equemene | SFMIN = DLAMCH( 'S' ) |
| 292 | 1 | equemene | SMLNUM = SFMIN / EPS |
| 293 | 1 | equemene | BIGNUM = ONE / SMLNUM |
| 294 | 1 | equemene | CALL DLABAD( SMLNUM, BIGNUM ) |
| 295 | 1 | equemene | * |
| 296 | 1 | equemene | * Scale A if max entry outside range [SMLNUM,BIGNUM]. |
| 297 | 1 | equemene | * |
| 298 | 1 | equemene | ANRM = DLANGE( 'M', M, N, A, LDA, WORK ) |
| 299 | 1 | equemene | IASCL = 0 |
| 300 | 1 | equemene | IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN |
| 301 | 1 | equemene | * |
| 302 | 1 | equemene | * Scale matrix norm up to SMLNUM. |
| 303 | 1 | equemene | * |
| 304 | 1 | equemene | CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) |
| 305 | 1 | equemene | IASCL = 1 |
| 306 | 1 | equemene | ELSE IF( ANRM.GT.BIGNUM ) THEN |
| 307 | 1 | equemene | * |
| 308 | 1 | equemene | * Scale matrix norm down to BIGNUM. |
| 309 | 1 | equemene | * |
| 310 | 1 | equemene | CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) |
| 311 | 1 | equemene | IASCL = 2 |
| 312 | 1 | equemene | ELSE IF( ANRM.EQ.ZERO ) THEN |
| 313 | 1 | equemene | * |
| 314 | 1 | equemene | * Matrix all zero. Return zero solution. |
| 315 | 1 | equemene | * |
| 316 | 1 | equemene | CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) |
| 317 | 1 | equemene | CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 ) |
| 318 | 1 | equemene | RANK = 0 |
| 319 | 1 | equemene | GO TO 10 |
| 320 | 1 | equemene | END IF |
| 321 | 1 | equemene | * |
| 322 | 1 | equemene | * Scale B if max entry outside range [SMLNUM,BIGNUM]. |
| 323 | 1 | equemene | * |
| 324 | 1 | equemene | BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK ) |
| 325 | 1 | equemene | IBSCL = 0 |
| 326 | 1 | equemene | IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN |
| 327 | 1 | equemene | * |
| 328 | 1 | equemene | * Scale matrix norm up to SMLNUM. |
| 329 | 1 | equemene | * |
| 330 | 1 | equemene | CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) |
| 331 | 1 | equemene | IBSCL = 1 |
| 332 | 1 | equemene | ELSE IF( BNRM.GT.BIGNUM ) THEN |
| 333 | 1 | equemene | * |
| 334 | 1 | equemene | * Scale matrix norm down to BIGNUM. |
| 335 | 1 | equemene | * |
| 336 | 1 | equemene | CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) |
| 337 | 1 | equemene | IBSCL = 2 |
| 338 | 1 | equemene | END IF |
| 339 | 1 | equemene | * |
| 340 | 1 | equemene | * If M < N make sure certain entries of B are zero. |
| 341 | 1 | equemene | * |
| 342 | 1 | equemene | IF( M.LT.N ) |
| 343 | 1 | equemene | $ CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) |
| 344 | 1 | equemene | * |
| 345 | 1 | equemene | * Overdetermined case. |
| 346 | 1 | equemene | * |
| 347 | 1 | equemene | IF( M.GE.N ) THEN |
| 348 | 1 | equemene | * |
| 349 | 1 | equemene | * Path 1 - overdetermined or exactly determined. |
| 350 | 1 | equemene | * |
| 351 | 1 | equemene | MM = M |
| 352 | 1 | equemene | IF( M.GE.MNTHR ) THEN |
| 353 | 1 | equemene | * |
| 354 | 1 | equemene | * Path 1a - overdetermined, with many more rows than columns. |
| 355 | 1 | equemene | * |
| 356 | 1 | equemene | MM = N |
| 357 | 1 | equemene | ITAU = 1 |
| 358 | 1 | equemene | NWORK = ITAU + N |
| 359 | 1 | equemene | * |
| 360 | 1 | equemene | * Compute A=Q*R. |
| 361 | 1 | equemene | * (Workspace: need 2*N, prefer N+N*NB) |
| 362 | 1 | equemene | * |
| 363 | 1 | equemene | CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
| 364 | 1 | equemene | $ LWORK-NWORK+1, INFO ) |
| 365 | 1 | equemene | * |
| 366 | 1 | equemene | * Multiply B by transpose(Q). |
| 367 | 1 | equemene | * (Workspace: need N+NRHS, prefer N+NRHS*NB) |
| 368 | 1 | equemene | * |
| 369 | 1 | equemene | CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B, |
| 370 | 1 | equemene | $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) |
| 371 | 1 | equemene | * |
| 372 | 1 | equemene | * Zero out below R. |
| 373 | 1 | equemene | * |
| 374 | 1 | equemene | IF( N.GT.1 ) THEN |
| 375 | 1 | equemene | CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) |
| 376 | 1 | equemene | END IF |
| 377 | 1 | equemene | END IF |
| 378 | 1 | equemene | * |
| 379 | 1 | equemene | IE = 1 |
| 380 | 1 | equemene | ITAUQ = IE + N |
| 381 | 1 | equemene | ITAUP = ITAUQ + N |
| 382 | 1 | equemene | NWORK = ITAUP + N |
| 383 | 1 | equemene | * |
| 384 | 1 | equemene | * Bidiagonalize R in A. |
| 385 | 1 | equemene | * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB) |
| 386 | 1 | equemene | * |
| 387 | 1 | equemene | CALL DGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
| 388 | 1 | equemene | $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
| 389 | 1 | equemene | $ INFO ) |
| 390 | 1 | equemene | * |
| 391 | 1 | equemene | * Multiply B by transpose of left bidiagonalizing vectors of R. |
| 392 | 1 | equemene | * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB) |
| 393 | 1 | equemene | * |
| 394 | 1 | equemene | CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ), |
| 395 | 1 | equemene | $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) |
| 396 | 1 | equemene | * |
| 397 | 1 | equemene | * Solve the bidiagonal least squares problem. |
| 398 | 1 | equemene | * |
| 399 | 1 | equemene | CALL DLALSD( 'U', SMLSIZ, N, NRHS, S, WORK( IE ), B, LDB, |
| 400 | 1 | equemene | $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) |
| 401 | 1 | equemene | IF( INFO.NE.0 ) THEN |
| 402 | 1 | equemene | GO TO 10 |
| 403 | 1 | equemene | END IF |
| 404 | 1 | equemene | * |
| 405 | 1 | equemene | * Multiply B by right bidiagonalizing vectors of R. |
| 406 | 1 | equemene | * |
| 407 | 1 | equemene | CALL DORMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ), |
| 408 | 1 | equemene | $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) |
| 409 | 1 | equemene | * |
| 410 | 1 | equemene | ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+ |
| 411 | 1 | equemene | $ MAX( M, 2*M-4, NRHS, N-3*M, WLALSD ) ) THEN |
| 412 | 1 | equemene | * |
| 413 | 1 | equemene | * Path 2a - underdetermined, with many more columns than rows |
| 414 | 1 | equemene | * and sufficient workspace for an efficient algorithm. |
| 415 | 1 | equemene | * |
| 416 | 1 | equemene | LDWORK = M |
| 417 | 1 | equemene | IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ), |
| 418 | 1 | equemene | $ M*LDA+M+M*NRHS, 4*M+M*LDA+WLALSD ) )LDWORK = LDA |
| 419 | 1 | equemene | ITAU = 1 |
| 420 | 1 | equemene | NWORK = M + 1 |
| 421 | 1 | equemene | * |
| 422 | 1 | equemene | * Compute A=L*Q. |
| 423 | 1 | equemene | * (Workspace: need 2*M, prefer M+M*NB) |
| 424 | 1 | equemene | * |
| 425 | 1 | equemene | CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
| 426 | 1 | equemene | $ LWORK-NWORK+1, INFO ) |
| 427 | 1 | equemene | IL = NWORK |
| 428 | 1 | equemene | * |
| 429 | 1 | equemene | * Copy L to WORK(IL), zeroing out above its diagonal. |
| 430 | 1 | equemene | * |
| 431 | 1 | equemene | CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK ) |
| 432 | 1 | equemene | CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ), |
| 433 | 1 | equemene | $ LDWORK ) |
| 434 | 1 | equemene | IE = IL + LDWORK*M |
| 435 | 1 | equemene | ITAUQ = IE + M |
| 436 | 1 | equemene | ITAUP = ITAUQ + M |
| 437 | 1 | equemene | NWORK = ITAUP + M |
| 438 | 1 | equemene | * |
| 439 | 1 | equemene | * Bidiagonalize L in WORK(IL). |
| 440 | 1 | equemene | * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB) |
| 441 | 1 | equemene | * |
| 442 | 1 | equemene | CALL DGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ), |
| 443 | 1 | equemene | $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), |
| 444 | 1 | equemene | $ LWORK-NWORK+1, INFO ) |
| 445 | 1 | equemene | * |
| 446 | 1 | equemene | * Multiply B by transpose of left bidiagonalizing vectors of L. |
| 447 | 1 | equemene | * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) |
| 448 | 1 | equemene | * |
| 449 | 1 | equemene | CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK, |
| 450 | 1 | equemene | $ WORK( ITAUQ ), B, LDB, WORK( NWORK ), |
| 451 | 1 | equemene | $ LWORK-NWORK+1, INFO ) |
| 452 | 1 | equemene | * |
| 453 | 1 | equemene | * Solve the bidiagonal least squares problem. |
| 454 | 1 | equemene | * |
| 455 | 1 | equemene | CALL DLALSD( 'U', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB, |
| 456 | 1 | equemene | $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) |
| 457 | 1 | equemene | IF( INFO.NE.0 ) THEN |
| 458 | 1 | equemene | GO TO 10 |
| 459 | 1 | equemene | END IF |
| 460 | 1 | equemene | * |
| 461 | 1 | equemene | * Multiply B by right bidiagonalizing vectors of L. |
| 462 | 1 | equemene | * |
| 463 | 1 | equemene | CALL DORMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK, |
| 464 | 1 | equemene | $ WORK( ITAUP ), B, LDB, WORK( NWORK ), |
| 465 | 1 | equemene | $ LWORK-NWORK+1, INFO ) |
| 466 | 1 | equemene | * |
| 467 | 1 | equemene | * Zero out below first M rows of B. |
| 468 | 1 | equemene | * |
| 469 | 1 | equemene | CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB ) |
| 470 | 1 | equemene | NWORK = ITAU + M |
| 471 | 1 | equemene | * |
| 472 | 1 | equemene | * Multiply transpose(Q) by B. |
| 473 | 1 | equemene | * (Workspace: need M+NRHS, prefer M+NRHS*NB) |
| 474 | 1 | equemene | * |
| 475 | 1 | equemene | CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B, |
| 476 | 1 | equemene | $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) |
| 477 | 1 | equemene | * |
| 478 | 1 | equemene | ELSE |
| 479 | 1 | equemene | * |
| 480 | 1 | equemene | * Path 2 - remaining underdetermined cases. |
| 481 | 1 | equemene | * |
| 482 | 1 | equemene | IE = 1 |
| 483 | 1 | equemene | ITAUQ = IE + M |
| 484 | 1 | equemene | ITAUP = ITAUQ + M |
| 485 | 1 | equemene | NWORK = ITAUP + M |
| 486 | 1 | equemene | * |
| 487 | 1 | equemene | * Bidiagonalize A. |
| 488 | 1 | equemene | * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) |
| 489 | 1 | equemene | * |
| 490 | 1 | equemene | CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
| 491 | 1 | equemene | $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
| 492 | 1 | equemene | $ INFO ) |
| 493 | 1 | equemene | * |
| 494 | 1 | equemene | * Multiply B by transpose of left bidiagonalizing vectors. |
| 495 | 1 | equemene | * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB) |
| 496 | 1 | equemene | * |
| 497 | 1 | equemene | CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ), |
| 498 | 1 | equemene | $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) |
| 499 | 1 | equemene | * |
| 500 | 1 | equemene | * Solve the bidiagonal least squares problem. |
| 501 | 1 | equemene | * |
| 502 | 1 | equemene | CALL DLALSD( 'L', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB, |
| 503 | 1 | equemene | $ RCOND, RANK, WORK( NWORK ), IWORK, INFO ) |
| 504 | 1 | equemene | IF( INFO.NE.0 ) THEN |
| 505 | 1 | equemene | GO TO 10 |
| 506 | 1 | equemene | END IF |
| 507 | 1 | equemene | * |
| 508 | 1 | equemene | * Multiply B by right bidiagonalizing vectors of A. |
| 509 | 1 | equemene | * |
| 510 | 1 | equemene | CALL DORMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ), |
| 511 | 1 | equemene | $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) |
| 512 | 1 | equemene | * |
| 513 | 1 | equemene | END IF |
| 514 | 1 | equemene | * |
| 515 | 1 | equemene | * Undo scaling. |
| 516 | 1 | equemene | * |
| 517 | 1 | equemene | IF( IASCL.EQ.1 ) THEN |
| 518 | 1 | equemene | CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) |
| 519 | 1 | equemene | CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, |
| 520 | 1 | equemene | $ INFO ) |
| 521 | 1 | equemene | ELSE IF( IASCL.EQ.2 ) THEN |
| 522 | 1 | equemene | CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) |
| 523 | 1 | equemene | CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, |
| 524 | 1 | equemene | $ INFO ) |
| 525 | 1 | equemene | END IF |
| 526 | 1 | equemene | IF( IBSCL.EQ.1 ) THEN |
| 527 | 1 | equemene | CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) |
| 528 | 1 | equemene | ELSE IF( IBSCL.EQ.2 ) THEN |
| 529 | 1 | equemene | CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) |
| 530 | 1 | equemene | END IF |
| 531 | 1 | equemene | * |
| 532 | 1 | equemene | 10 CONTINUE |
| 533 | 1 | equemene | WORK( 1 ) = MAXWRK |
| 534 | 1 | equemene | IWORK( 1 ) = LIWORK |
| 535 | 1 | equemene | RETURN |
| 536 | 1 | equemene | * |
| 537 | 1 | equemene | * End of DGELSD |
| 538 | 1 | equemene | * |
| 539 | 1 | equemene | END |