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