LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
sgelsy.f
Go to the documentation of this file.
1 *> \brief <b> SGELSY 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 SGELSY + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgelsy.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgelsy.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgelsy.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SGELSY( 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 * REAL RCOND
27 * ..
28 * .. Array Arguments ..
29 * INTEGER JPVT( * )
30 * REAL A( LDA, * ), B( LDB, * ), WORK( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SGELSY 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 REAL 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 REAL 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 REAL
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 REAL 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 SGEQP3, STZRZF, STZRQF, SORMQR,
169 *> and SORMRZ.
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 *> \date November 2011
193 *
194 *> \ingroup realGEsolve
195 *
196 *> \par Contributors:
197 * ==================
198 *>
199 *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n
200 *> E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
201 *> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
202 *>
203 * =====================================================================
204  SUBROUTINE sgelsy( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
205  $ work, lwork, info )
206 *
207 * -- LAPACK driver routine (version 3.4.0) --
208 * -- LAPACK is a software package provided by Univ. of Tennessee, --
209 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
210 * November 2011
211 *
212 * .. Scalar Arguments ..
213  INTEGER info, lda, ldb, lwork, m, n, nrhs, rank
214  REAL rcond
215 * ..
216 * .. Array Arguments ..
217  INTEGER jpvt( * )
218  REAL a( lda, * ), b( ldb, * ), work( * )
219 * ..
220 *
221 * =====================================================================
222 *
223 * .. Parameters ..
224  INTEGER imax, imin
225  parameter( imax = 1, imin = 2 )
226  REAL zero, one
227  parameter( zero = 0.0e+0, one = 1.0e+0 )
228 * ..
229 * .. Local Scalars ..
230  LOGICAL lquery
231  INTEGER i, iascl, ibscl, ismax, ismin, j, lwkmin,
232  $ lwkopt, mn, nb, nb1, nb2, nb3, nb4
233  REAL anrm, bignum, bnrm, c1, c2, s1, s2, smax,
234  $ smaxpr, smin, sminpr, smlnum, wsize
235 * ..
236 * .. External Functions ..
237  INTEGER ilaenv
238  REAL slamch, slange
239  EXTERNAL ilaenv, slamch, slange
240 * ..
241 * .. External Subroutines ..
242  EXTERNAL scopy, sgeqp3, slabad, slaic1, slascl, slaset,
244 * ..
245 * .. Intrinsic Functions ..
246  INTRINSIC abs, max, min
247 * ..
248 * .. Executable Statements ..
249 *
250  mn = min( m, n )
251  ismin = mn + 1
252  ismax = 2*mn + 1
253 *
254 * Test the input arguments.
255 *
256  info = 0
257  lquery = ( lwork.EQ.-1 )
258  IF( m.LT.0 ) THEN
259  info = -1
260  ELSE IF( n.LT.0 ) THEN
261  info = -2
262  ELSE IF( nrhs.LT.0 ) THEN
263  info = -3
264  ELSE IF( lda.LT.max( 1, m ) ) THEN
265  info = -5
266  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
267  info = -7
268  END IF
269 *
270 * Figure out optimal block size
271 *
272  IF( info.EQ.0 ) THEN
273  IF( mn.EQ.0 .OR. nrhs.EQ.0 ) THEN
274  lwkmin = 1
275  lwkopt = 1
276  ELSE
277  nb1 = ilaenv( 1, 'SGEQRF', ' ', m, n, -1, -1 )
278  nb2 = ilaenv( 1, 'SGERQF', ' ', m, n, -1, -1 )
279  nb3 = ilaenv( 1, 'SORMQR', ' ', m, n, nrhs, -1 )
280  nb4 = ilaenv( 1, 'SORMRQ', ' ', m, n, nrhs, -1 )
281  nb = max( nb1, nb2, nb3, nb4 )
282  lwkmin = mn + max( 2*mn, n + 1, mn + nrhs )
283  lwkopt = max( lwkmin,
284  $ mn + 2*n + nb*( n + 1 ), 2*mn + nb*nrhs )
285  END IF
286  work( 1 ) = lwkopt
287 *
288  IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
289  info = -12
290  END IF
291  END IF
292 *
293  IF( info.NE.0 ) THEN
294  CALL xerbla( 'SGELSY', -info )
295  return
296  ELSE IF( lquery ) THEN
297  return
298  END IF
299 *
300 * Quick return if possible
301 *
302  IF( mn.EQ.0 .OR. nrhs.EQ.0 ) THEN
303  rank = 0
304  return
305  END IF
306 *
307 * Get machine parameters
308 *
309  smlnum = slamch( 'S' ) / slamch( 'P' )
310  bignum = one / smlnum
311  CALL slabad( smlnum, bignum )
312 *
313 * Scale A, B if max entries outside range [SMLNUM,BIGNUM]
314 *
315  anrm = slange( 'M', m, n, a, lda, work )
316  iascl = 0
317  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
318 *
319 * Scale matrix norm up to SMLNUM
320 *
321  CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
322  iascl = 1
323  ELSE IF( anrm.GT.bignum ) THEN
324 *
325 * Scale matrix norm down to BIGNUM
326 *
327  CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
328  iascl = 2
329  ELSE IF( anrm.EQ.zero ) THEN
330 *
331 * Matrix all zero. Return zero solution.
332 *
333  CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
334  rank = 0
335  go to 70
336  END IF
337 *
338  bnrm = slange( 'M', m, nrhs, b, ldb, work )
339  ibscl = 0
340  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
341 *
342 * Scale matrix norm up to SMLNUM
343 *
344  CALL slascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
345  ibscl = 1
346  ELSE IF( bnrm.GT.bignum ) THEN
347 *
348 * Scale matrix norm down to BIGNUM
349 *
350  CALL slascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
351  ibscl = 2
352  END IF
353 *
354 * Compute QR factorization with column pivoting of A:
355 * A * P = Q * R
356 *
357  CALL sgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
358  $ lwork-mn, info )
359  wsize = mn + work( mn+1 )
360 *
361 * workspace: MN+2*N+NB*(N+1).
362 * Details of Householder rotations stored in WORK(1:MN).
363 *
364 * Determine RANK using incremental condition estimation
365 *
366  work( ismin ) = one
367  work( ismax ) = one
368  smax = abs( a( 1, 1 ) )
369  smin = smax
370  IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
371  rank = 0
372  CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
373  go to 70
374  ELSE
375  rank = 1
376  END IF
377 *
378  10 continue
379  IF( rank.LT.mn ) THEN
380  i = rank + 1
381  CALL slaic1( imin, rank, work( ismin ), smin, a( 1, i ),
382  $ a( i, i ), sminpr, s1, c1 )
383  CALL slaic1( imax, rank, work( ismax ), smax, a( 1, i ),
384  $ a( i, i ), smaxpr, s2, c2 )
385 *
386  IF( smaxpr*rcond.LE.sminpr ) THEN
387  DO 20 i = 1, rank
388  work( ismin+i-1 ) = s1*work( ismin+i-1 )
389  work( ismax+i-1 ) = s2*work( ismax+i-1 )
390  20 continue
391  work( ismin+rank ) = c1
392  work( ismax+rank ) = c2
393  smin = sminpr
394  smax = smaxpr
395  rank = rank + 1
396  go to 10
397  END IF
398  END IF
399 *
400 * workspace: 3*MN.
401 *
402 * Logically partition R = [ R11 R12 ]
403 * [ 0 R22 ]
404 * where R11 = R(1:RANK,1:RANK)
405 *
406 * [R11,R12] = [ T11, 0 ] * Y
407 *
408  IF( rank.LT.n )
409  $ CALL stzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
410  $ lwork-2*mn, info )
411 *
412 * workspace: 2*MN.
413 * Details of Householder rotations stored in WORK(MN+1:2*MN)
414 *
415 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
416 *
417  CALL sormqr( 'Left', 'Transpose', m, nrhs, mn, a, lda, work( 1 ),
418  $ b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
419  wsize = max( wsize, 2*mn+work( 2*mn+1 ) )
420 *
421 * workspace: 2*MN+NB*NRHS.
422 *
423 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
424 *
425  CALL strsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
426  $ nrhs, one, a, lda, b, ldb )
427 *
428  DO 40 j = 1, nrhs
429  DO 30 i = rank + 1, n
430  b( i, j ) = zero
431  30 continue
432  40 continue
433 *
434 * B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
435 *
436  IF( rank.LT.n ) THEN
437  CALL sormrz( 'Left', 'Transpose', n, nrhs, rank, n-rank, a,
438  $ lda, work( mn+1 ), b, ldb, work( 2*mn+1 ),
439  $ lwork-2*mn, info )
440  END IF
441 *
442 * workspace: 2*MN+NRHS.
443 *
444 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
445 *
446  DO 60 j = 1, nrhs
447  DO 50 i = 1, n
448  work( jpvt( i ) ) = b( i, j )
449  50 continue
450  CALL scopy( n, work( 1 ), 1, b( 1, j ), 1 )
451  60 continue
452 *
453 * workspace: N.
454 *
455 * Undo scaling
456 *
457  IF( iascl.EQ.1 ) THEN
458  CALL slascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
459  CALL slascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
460  $ info )
461  ELSE IF( iascl.EQ.2 ) THEN
462  CALL slascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
463  CALL slascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
464  $ info )
465  END IF
466  IF( ibscl.EQ.1 ) THEN
467  CALL slascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
468  ELSE IF( ibscl.EQ.2 ) THEN
469  CALL slascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
470  END IF
471 *
472  70 continue
473  work( 1 ) = lwkopt
474 *
475  return
476 *
477 * End of SGELSY
478 *
479  END