1: *> \brief <b> DGELSY solves overdetermined or underdetermined systems for GE matrices</b>
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DGELSY + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelsy.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelsy.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelsy.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
22: * WORK, LWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
26: * DOUBLE PRECISION RCOND
27: * ..
28: * .. Array Arguments ..
29: * INTEGER JPVT( * )
30: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DGELSY computes the minimum-norm solution to a real linear least
40: *> squares problem:
41: *> minimize || A * X - B ||
42: *> using a complete orthogonal factorization of A. A is an M-by-N
43: *> matrix which may be rank-deficient.
44: *>
45: *> Several right hand side vectors b and solution vectors x can be
46: *> handled in a single call; they are stored as the columns of the
47: *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
48: *> matrix X.
49: *>
50: *> The routine first computes a QR factorization with column pivoting:
51: *> A * P = Q * [ R11 R12 ]
52: *> [ 0 R22 ]
53: *> with R11 defined as the largest leading submatrix whose estimated
54: *> condition number is less than 1/RCOND. The order of R11, RANK,
55: *> is the effective rank of A.
56: *>
57: *> Then, R22 is considered to be negligible, and R12 is annihilated
58: *> by orthogonal transformations from the right, arriving at the
59: *> complete orthogonal factorization:
60: *> A * P = Q * [ T11 0 ] * Z
61: *> [ 0 0 ]
62: *> The minimum-norm solution is then
63: *> X = P * Z**T [ inv(T11)*Q1**T*B ]
64: *> [ 0 ]
65: *> where Q1 consists of the first RANK columns of Q.
66: *>
67: *> This routine is basically identical to the original xGELSX except
68: *> three differences:
69: *> o The call to the subroutine xGEQPF has been substituted by the
70: *> the call to the subroutine xGEQP3. This subroutine is a Blas-3
71: *> version of the QR factorization with column pivoting.
72: *> o Matrix B (the right hand side) is updated with Blas-3.
73: *> o The permutation of matrix B (the right hand side) is faster and
74: *> more simple.
75: *> \endverbatim
76: *
77: * Arguments:
78: * ==========
79: *
80: *> \param[in] M
81: *> \verbatim
82: *> M is INTEGER
83: *> The number of rows of the matrix A. M >= 0.
84: *> \endverbatim
85: *>
86: *> \param[in] N
87: *> \verbatim
88: *> N is INTEGER
89: *> The number of columns of the matrix A. N >= 0.
90: *> \endverbatim
91: *>
92: *> \param[in] NRHS
93: *> \verbatim
94: *> NRHS is INTEGER
95: *> The number of right hand sides, i.e., the number of
96: *> columns of matrices B and X. NRHS >= 0.
97: *> \endverbatim
98: *>
99: *> \param[in,out] A
100: *> \verbatim
101: *> A is DOUBLE PRECISION array, dimension (LDA,N)
102: *> On entry, the M-by-N matrix A.
103: *> On exit, A has been overwritten by details of its
104: *> complete orthogonal factorization.
105: *> \endverbatim
106: *>
107: *> \param[in] LDA
108: *> \verbatim
109: *> LDA is INTEGER
110: *> The leading dimension of the array A. LDA >= max(1,M).
111: *> \endverbatim
112: *>
113: *> \param[in,out] B
114: *> \verbatim
115: *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
116: *> On entry, the M-by-NRHS right hand side matrix B.
117: *> On exit, the N-by-NRHS solution matrix X.
118: *> \endverbatim
119: *>
120: *> \param[in] LDB
121: *> \verbatim
122: *> LDB is INTEGER
123: *> The leading dimension of the array B. LDB >= max(1,M,N).
124: *> \endverbatim
125: *>
126: *> \param[in,out] JPVT
127: *> \verbatim
128: *> JPVT is INTEGER array, dimension (N)
129: *> On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
130: *> to the front of AP, otherwise column i is a free column.
131: *> On exit, if JPVT(i) = k, then the i-th column of AP
132: *> was the k-th column of A.
133: *> \endverbatim
134: *>
135: *> \param[in] RCOND
136: *> \verbatim
137: *> RCOND is DOUBLE PRECISION
138: *> RCOND is used to determine the effective rank of A, which
139: *> is defined as the order of the largest leading triangular
140: *> submatrix R11 in the QR factorization with pivoting of A,
141: *> whose estimated condition number < 1/RCOND.
142: *> \endverbatim
143: *>
144: *> \param[out] RANK
145: *> \verbatim
146: *> RANK is INTEGER
147: *> The effective rank of A, i.e., the order of the submatrix
148: *> R11. This is the same as the order of the submatrix T11
149: *> in the complete orthogonal factorization of A.
150: *> \endverbatim
151: *>
152: *> \param[out] WORK
153: *> \verbatim
154: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
155: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
156: *> \endverbatim
157: *>
158: *> \param[in] LWORK
159: *> \verbatim
160: *> LWORK is INTEGER
161: *> The dimension of the array WORK.
162: *> The unblocked strategy requires that:
163: *> LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ),
164: *> where MN = min( M, N ).
165: *> The block algorithm requires that:
166: *> LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),
167: *> where NB is an upper bound on the blocksize returned
168: *> by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,
169: *> and DORMRZ.
170: *>
171: *> If LWORK = -1, then a workspace query is assumed; the routine
172: *> only calculates the optimal size of the WORK array, returns
173: *> this value as the first entry of the WORK array, and no error
174: *> message related to LWORK is issued by XERBLA.
175: *> \endverbatim
176: *>
177: *> \param[out] INFO
178: *> \verbatim
179: *> INFO is INTEGER
180: *> = 0: successful exit
181: *> < 0: If INFO = -i, the i-th argument had an illegal value.
182: *> \endverbatim
183: *
184: * Authors:
185: * ========
186: *
187: *> \author Univ. of Tennessee
188: *> \author Univ. of California Berkeley
189: *> \author Univ. of Colorado Denver
190: *> \author NAG Ltd.
191: *
192: *> \ingroup doubleGEsolve
193: *
194: *> \par Contributors:
195: * ==================
196: *>
197: *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n
198: *> E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
199: *> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
200: *>
201: * =====================================================================
202: SUBROUTINE DGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
203: $ WORK, LWORK, INFO )
204: *
205: * -- LAPACK driver routine --
206: * -- LAPACK is a software package provided by Univ. of Tennessee, --
207: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
208: *
209: * .. Scalar Arguments ..
210: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
211: DOUBLE PRECISION RCOND
212: * ..
213: * .. Array Arguments ..
214: INTEGER JPVT( * )
215: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
216: * ..
217: *
218: * =====================================================================
219: *
220: * .. Parameters ..
221: INTEGER IMAX, IMIN
222: PARAMETER ( IMAX = 1, IMIN = 2 )
223: DOUBLE PRECISION ZERO, ONE
224: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
225: * ..
226: * .. Local Scalars ..
227: LOGICAL LQUERY
228: INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKMIN,
229: $ LWKOPT, MN, NB, NB1, NB2, NB3, NB4
230: DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
231: $ SMAXPR, SMIN, SMINPR, SMLNUM, WSIZE
232: * ..
233: * .. External Functions ..
234: INTEGER ILAENV
235: DOUBLE PRECISION DLAMCH, DLANGE
236: EXTERNAL ILAENV, DLAMCH, DLANGE
237: * ..
238: * .. External Subroutines ..
239: EXTERNAL DCOPY, DGEQP3, DLABAD, DLAIC1, DLASCL, DLASET,
240: $ DORMQR, DORMRZ, DTRSM, DTZRZF, XERBLA
241: * ..
242: * .. Intrinsic Functions ..
243: INTRINSIC ABS, MAX, MIN
244: * ..
245: * .. Executable Statements ..
246: *
247: MN = MIN( M, N )
248: ISMIN = MN + 1
249: ISMAX = 2*MN + 1
250: *
251: * Test the input arguments.
252: *
253: INFO = 0
254: LQUERY = ( LWORK.EQ.-1 )
255: IF( M.LT.0 ) THEN
256: INFO = -1
257: ELSE IF( N.LT.0 ) THEN
258: INFO = -2
259: ELSE IF( NRHS.LT.0 ) THEN
260: INFO = -3
261: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
262: INFO = -5
263: ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
264: INFO = -7
265: END IF
266: *
267: * Figure out optimal block size
268: *
269: IF( INFO.EQ.0 ) THEN
270: IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN
271: LWKMIN = 1
272: LWKOPT = 1
273: ELSE
274: NB1 = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
275: NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 )
276: NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, NRHS, -1 )
277: NB4 = ILAENV( 1, 'DORMRQ', ' ', M, N, NRHS, -1 )
278: NB = MAX( NB1, NB2, NB3, NB4 )
279: LWKMIN = MN + MAX( 2*MN, N + 1, MN + NRHS )
280: LWKOPT = MAX( LWKMIN,
281: $ MN + 2*N + NB*( N + 1 ), 2*MN + NB*NRHS )
282: END IF
283: WORK( 1 ) = LWKOPT
284: *
285: IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
286: INFO = -12
287: END IF
288: END IF
289: *
290: IF( INFO.NE.0 ) THEN
291: CALL XERBLA( 'DGELSY', -INFO )
292: RETURN
293: ELSE IF( LQUERY ) THEN
294: RETURN
295: END IF
296: *
297: * Quick return if possible
298: *
299: IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN
300: RANK = 0
301: RETURN
302: END IF
303: *
304: * Get machine parameters
305: *
306: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
307: BIGNUM = ONE / SMLNUM
308: CALL DLABAD( SMLNUM, BIGNUM )
309: *
310: * Scale A, B if max entries outside range [SMLNUM,BIGNUM]
311: *
312: ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
313: IASCL = 0
314: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
315: *
316: * Scale matrix norm up to SMLNUM
317: *
318: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
319: IASCL = 1
320: ELSE IF( ANRM.GT.BIGNUM ) THEN
321: *
322: * Scale matrix norm down to BIGNUM
323: *
324: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
325: IASCL = 2
326: ELSE IF( ANRM.EQ.ZERO ) THEN
327: *
328: * Matrix all zero. Return zero solution.
329: *
330: CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
331: RANK = 0
332: GO TO 70
333: END IF
334: *
335: BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
336: IBSCL = 0
337: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
338: *
339: * Scale matrix norm up to SMLNUM
340: *
341: CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
342: IBSCL = 1
343: ELSE IF( BNRM.GT.BIGNUM ) THEN
344: *
345: * Scale matrix norm down to BIGNUM
346: *
347: CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
348: IBSCL = 2
349: END IF
350: *
351: * Compute QR factorization with column pivoting of A:
352: * A * P = Q * R
353: *
354: CALL DGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ),
355: $ LWORK-MN, INFO )
356: WSIZE = MN + WORK( MN+1 )
357: *
358: * workspace: MN+2*N+NB*(N+1).
359: * Details of Householder rotations stored in WORK(1:MN).
360: *
361: * Determine RANK using incremental condition estimation
362: *
363: WORK( ISMIN ) = ONE
364: WORK( ISMAX ) = ONE
365: SMAX = ABS( A( 1, 1 ) )
366: SMIN = SMAX
367: IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
368: RANK = 0
369: CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
370: GO TO 70
371: ELSE
372: RANK = 1
373: END IF
374: *
375: 10 CONTINUE
376: IF( RANK.LT.MN ) THEN
377: I = RANK + 1
378: CALL DLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
379: $ A( I, I ), SMINPR, S1, C1 )
380: CALL DLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
381: $ A( I, I ), SMAXPR, S2, C2 )
382: *
383: IF( SMAXPR*RCOND.LE.SMINPR ) THEN
384: DO 20 I = 1, RANK
385: WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
386: WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
387: 20 CONTINUE
388: WORK( ISMIN+RANK ) = C1
389: WORK( ISMAX+RANK ) = C2
390: SMIN = SMINPR
391: SMAX = SMAXPR
392: RANK = RANK + 1
393: GO TO 10
394: END IF
395: END IF
396: *
397: * workspace: 3*MN.
398: *
399: * Logically partition R = [ R11 R12 ]
400: * [ 0 R22 ]
401: * where R11 = R(1:RANK,1:RANK)
402: *
403: * [R11,R12] = [ T11, 0 ] * Y
404: *
405: IF( RANK.LT.N )
406: $ CALL DTZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ),
407: $ LWORK-2*MN, INFO )
408: *
409: * workspace: 2*MN.
410: * Details of Householder rotations stored in WORK(MN+1:2*MN)
411: *
412: * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
413: *
414: CALL DORMQR( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ),
415: $ B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO )
416: WSIZE = MAX( WSIZE, 2*MN+WORK( 2*MN+1 ) )
417: *
418: * workspace: 2*MN+NB*NRHS.
419: *
420: * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
421: *
422: CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
423: $ NRHS, ONE, A, LDA, B, LDB )
424: *
425: DO 40 J = 1, NRHS
426: DO 30 I = RANK + 1, N
427: B( I, J ) = ZERO
428: 30 CONTINUE
429: 40 CONTINUE
430: *
431: * B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
432: *
433: IF( RANK.LT.N ) THEN
434: CALL DORMRZ( 'Left', 'Transpose', N, NRHS, RANK, N-RANK, A,
435: $ LDA, WORK( MN+1 ), B, LDB, WORK( 2*MN+1 ),
436: $ LWORK-2*MN, INFO )
437: END IF
438: *
439: * workspace: 2*MN+NRHS.
440: *
441: * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
442: *
443: DO 60 J = 1, NRHS
444: DO 50 I = 1, N
445: WORK( JPVT( I ) ) = B( I, J )
446: 50 CONTINUE
447: CALL DCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 )
448: 60 CONTINUE
449: *
450: * workspace: N.
451: *
452: * Undo scaling
453: *
454: IF( IASCL.EQ.1 ) THEN
455: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
456: CALL DLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
457: $ INFO )
458: ELSE IF( IASCL.EQ.2 ) THEN
459: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
460: CALL DLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
461: $ INFO )
462: END IF
463: IF( IBSCL.EQ.1 ) THEN
464: CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
465: ELSE IF( IBSCL.EQ.2 ) THEN
466: CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
467: END IF
468: *
469: 70 CONTINUE
470: WORK( 1 ) = LWKOPT
471: *
472: RETURN
473: *
474: * End of DGELSY
475: *
476: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgelsy.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack