1 /* ./src_f77/sgelsx.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__0 = 0;
11 static real c_b13 = 0.f;
12 static integer c__2 = 2;
13 static integer c__1 = 1;
14 static real c_b36 = 1.f;
15
sgelsx_(integer * m,integer * n,integer * nrhs,real * a,integer * lda,real * b,integer * ldb,integer * jpvt,real * rcond,integer * rank,real * work,integer * info)16 /* Subroutine */ int sgelsx_(integer *m, integer *n, integer *nrhs, real *a,
17 integer *lda, real *b, integer *ldb, integer *jpvt, real *rcond,
18 integer *rank, real *work, integer *info)
19 {
20 /* System generated locals */
21 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
22 real r__1;
23
24 /* Local variables */
25 static integer i__, j, k;
26 static real c1, c2, s1, s2, t1, t2;
27 static integer mn;
28 static real anrm, bnrm, smin, smax;
29 static integer iascl, ibscl, ismin, ismax;
30 extern /* Subroutine */ int strsm_(char *, char *, char *, char *,
31 integer *, integer *, real *, real *, integer *, real *, integer *
32 , ftnlen, ftnlen, ftnlen, ftnlen), slaic1_(integer *, integer *,
33 real *, real *, real *, real *, real *, real *, real *), sorm2r_(
34 char *, char *, integer *, integer *, integer *, real *, integer *
35 , real *, real *, integer *, real *, integer *, ftnlen, ftnlen),
36 slabad_(real *, real *);
37 extern doublereal slamch_(char *, ftnlen), slange_(char *, integer *,
38 integer *, real *, integer *, real *, ftnlen);
39 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
40 static real bignum;
41 extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *,
42 real *, integer *, integer *, real *, integer *, integer *,
43 ftnlen), sgeqpf_(integer *, integer *, real *, integer *, integer
44 *, real *, real *, integer *), slaset_(char *, integer *, integer
45 *, real *, real *, real *, integer *, ftnlen);
46 static real sminpr, smaxpr, smlnum;
47 extern /* Subroutine */ int slatzm_(char *, integer *, integer *, real *,
48 integer *, real *, real *, real *, integer *, real *, ftnlen),
49 stzrqf_(integer *, integer *, real *, integer *, real *, integer *
50 );
51
52
53 /* -- LAPACK driver routine (version 3.0) -- */
54 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
55 /* Courant Institute, Argonne National Lab, and Rice University */
56 /* March 31, 1993 */
57
58 /* .. Scalar Arguments .. */
59 /* .. */
60 /* .. Array Arguments .. */
61 /* .. */
62
63 /* Purpose */
64 /* ======= */
65
66 /* This routine is deprecated and has been replaced by routine SGELSY. */
67
68 /* SGELSX computes the minimum-norm solution to a real linear least */
69 /* squares problem: */
70 /* minimize || A * X - B || */
71 /* using a complete orthogonal factorization of A. A is an M-by-N */
72 /* matrix which may be rank-deficient. */
73
74 /* Several right hand side vectors b and solution vectors x can be */
75 /* handled in a single call; they are stored as the columns of the */
76 /* M-by-NRHS right hand side matrix B and the N-by-NRHS solution */
77 /* matrix X. */
78
79 /* The routine first computes a QR factorization with column pivoting: */
80 /* A * P = Q * [ R11 R12 ] */
81 /* [ 0 R22 ] */
82 /* with R11 defined as the largest leading submatrix whose estimated */
83 /* condition number is less than 1/RCOND. The order of R11, RANK, */
84 /* is the effective rank of A. */
85
86 /* Then, R22 is considered to be negligible, and R12 is annihilated */
87 /* by orthogonal transformations from the right, arriving at the */
88 /* complete orthogonal factorization: */
89 /* A * P = Q * [ T11 0 ] * Z */
90 /* [ 0 0 ] */
91 /* The minimum-norm solution is then */
92 /* X = P * Z' [ inv(T11)*Q1'*B ] */
93 /* [ 0 ] */
94 /* where Q1 consists of the first RANK columns of Q. */
95
96 /* Arguments */
97 /* ========= */
98
99 /* M (input) INTEGER */
100 /* The number of rows of the matrix A. M >= 0. */
101
102 /* N (input) INTEGER */
103 /* The number of columns of the matrix A. N >= 0. */
104
105 /* NRHS (input) INTEGER */
106 /* The number of right hand sides, i.e., the number of */
107 /* columns of matrices B and X. NRHS >= 0. */
108
109 /* A (input/output) REAL array, dimension (LDA,N) */
110 /* On entry, the M-by-N matrix A. */
111 /* On exit, A has been overwritten by details of its */
112 /* complete orthogonal factorization. */
113
114 /* LDA (input) INTEGER */
115 /* The leading dimension of the array A. LDA >= max(1,M). */
116
117 /* B (input/output) REAL array, dimension (LDB,NRHS) */
118 /* On entry, the M-by-NRHS right hand side matrix B. */
119 /* On exit, the N-by-NRHS solution matrix X. */
120 /* If m >= n and RANK = n, the residual sum-of-squares for */
121 /* the solution in the i-th column is given by the sum of */
122 /* squares of elements N+1:M in that column. */
123
124 /* LDB (input) INTEGER */
125 /* The leading dimension of the array B. LDB >= max(1,M,N). */
126
127 /* JPVT (input/output) INTEGER array, dimension (N) */
128 /* On entry, if JPVT(i) .ne. 0, the i-th column of A is an */
129 /* initial column, otherwise it is a free column. Before */
130 /* the QR factorization of A, all initial columns are */
131 /* permuted to the leading positions; only the remaining */
132 /* free columns are moved as a result of column pivoting */
133 /* during the factorization. */
134 /* On exit, if JPVT(i) = k, then the i-th column of A*P */
135 /* was the k-th column of A. */
136
137 /* RCOND (input) 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
143 /* RANK (output) INTEGER */
144 /* The effective rank of A, i.e., the order of the submatrix */
145 /* R11. This is the same as the order of the submatrix T11 */
146 /* in the complete orthogonal factorization of A. */
147
148 /* WORK (workspace) REAL array, dimension */
149 /* (max( min(M,N)+3*N, 2*min(M,N)+NRHS )), */
150
151 /* INFO (output) INTEGER */
152 /* = 0: successful exit */
153 /* < 0: if INFO = -i, the i-th argument had an illegal value */
154
155 /* ===================================================================== */
156
157 /* .. Parameters .. */
158 /* .. */
159 /* .. Local Scalars .. */
160 /* .. */
161 /* .. External Functions .. */
162 /* .. */
163 /* .. External Subroutines .. */
164 /* .. */
165 /* .. Intrinsic Functions .. */
166 /* .. */
167 /* .. Executable Statements .. */
168
169 /* Parameter adjustments */
170 a_dim1 = *lda;
171 a_offset = 1 + a_dim1;
172 a -= a_offset;
173 b_dim1 = *ldb;
174 b_offset = 1 + b_dim1;
175 b -= b_offset;
176 --jpvt;
177 --work;
178
179 /* Function Body */
180 mn = min(*m,*n);
181 ismin = mn + 1;
182 ismax = (mn << 1) + 1;
183
184 /* Test the input arguments. */
185
186 *info = 0;
187 if (*m < 0) {
188 *info = -1;
189 } else if (*n < 0) {
190 *info = -2;
191 } else if (*nrhs < 0) {
192 *info = -3;
193 } else if (*lda < max(1,*m)) {
194 *info = -5;
195 } else /* if(complicated condition) */ {
196 /* Computing MAX */
197 i__1 = max(1,*m);
198 if (*ldb < max(i__1,*n)) {
199 *info = -7;
200 }
201 }
202
203 if (*info != 0) {
204 i__1 = -(*info);
205 xerbla_("SGELSX", &i__1, (ftnlen)6);
206 return 0;
207 }
208
209 /* Quick return if possible */
210
211 /* Computing MIN */
212 i__1 = min(*m,*n);
213 if (min(i__1,*nrhs) == 0) {
214 *rank = 0;
215 return 0;
216 }
217
218 /* Get machine parameters */
219
220 smlnum = slamch_("S", (ftnlen)1) / slamch_("P", (ftnlen)1);
221 bignum = 1.f / smlnum;
222 slabad_(&smlnum, &bignum);
223
224 /* Scale A, B if max elements outside range [SMLNUM,BIGNUM] */
225
226 anrm = slange_("M", m, n, &a[a_offset], lda, &work[1], (ftnlen)1);
227 iascl = 0;
228 if (anrm > 0.f && anrm < smlnum) {
229
230 /* Scale matrix norm up to SMLNUM */
231
232 slascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda,
233 info, (ftnlen)1);
234 iascl = 1;
235 } else if (anrm > bignum) {
236
237 /* Scale matrix norm down to BIGNUM */
238
239 slascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda,
240 info, (ftnlen)1);
241 iascl = 2;
242 } else if (anrm == 0.f) {
243
244 /* Matrix all zero. Return zero solution. */
245
246 i__1 = max(*m,*n);
247 slaset_("F", &i__1, nrhs, &c_b13, &c_b13, &b[b_offset], ldb, (ftnlen)
248 1);
249 *rank = 0;
250 goto L100;
251 }
252
253 bnrm = slange_("M", m, nrhs, &b[b_offset], ldb, &work[1], (ftnlen)1);
254 ibscl = 0;
255 if (bnrm > 0.f && bnrm < smlnum) {
256
257 /* Scale matrix norm up to SMLNUM */
258
259 slascl_("G", &c__0, &c__0, &bnrm, &smlnum, m, nrhs, &b[b_offset], ldb,
260 info, (ftnlen)1);
261 ibscl = 1;
262 } else if (bnrm > bignum) {
263
264 /* Scale matrix norm down to BIGNUM */
265
266 slascl_("G", &c__0, &c__0, &bnrm, &bignum, m, nrhs, &b[b_offset], ldb,
267 info, (ftnlen)1);
268 ibscl = 2;
269 }
270
271 /* Compute QR factorization with column pivoting of A: */
272 /* A * P = Q * R */
273
274 sgeqpf_(m, n, &a[a_offset], lda, &jpvt[1], &work[1], &work[mn + 1], info);
275
276 /* workspace 3*N. Details of Householder rotations stored */
277 /* in WORK(1:MN). */
278
279 /* Determine RANK using incremental condition estimation */
280
281 work[ismin] = 1.f;
282 work[ismax] = 1.f;
283 smax = (r__1 = a[a_dim1 + 1], dabs(r__1));
284 smin = smax;
285 if ((r__1 = a[a_dim1 + 1], dabs(r__1)) == 0.f) {
286 *rank = 0;
287 i__1 = max(*m,*n);
288 slaset_("F", &i__1, nrhs, &c_b13, &c_b13, &b[b_offset], ldb, (ftnlen)
289 1);
290 goto L100;
291 } else {
292 *rank = 1;
293 }
294
295 L10:
296 if (*rank < mn) {
297 i__ = *rank + 1;
298 slaic1_(&c__2, rank, &work[ismin], &smin, &a[i__ * a_dim1 + 1], &a[
299 i__ + i__ * a_dim1], &sminpr, &s1, &c1);
300 slaic1_(&c__1, rank, &work[ismax], &smax, &a[i__ * a_dim1 + 1], &a[
301 i__ + i__ * a_dim1], &smaxpr, &s2, &c2);
302
303 if (smaxpr * *rcond <= sminpr) {
304 i__1 = *rank;
305 for (i__ = 1; i__ <= i__1; ++i__) {
306 work[ismin + i__ - 1] = s1 * work[ismin + i__ - 1];
307 work[ismax + i__ - 1] = s2 * work[ismax + i__ - 1];
308 /* L20: */
309 }
310 work[ismin + *rank] = c1;
311 work[ismax + *rank] = c2;
312 smin = sminpr;
313 smax = smaxpr;
314 ++(*rank);
315 goto L10;
316 }
317 }
318
319 /* Logically partition R = [ R11 R12 ] */
320 /* [ 0 R22 ] */
321 /* where R11 = R(1:RANK,1:RANK) */
322
323 /* [R11,R12] = [ T11, 0 ] * Y */
324
325 if (*rank < *n) {
326 stzrqf_(rank, n, &a[a_offset], lda, &work[mn + 1], info);
327 }
328
329 /* Details of Householder rotations stored in WORK(MN+1:2*MN) */
330
331 /* B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) */
332
333 sorm2r_("Left", "Transpose", m, nrhs, &mn, &a[a_offset], lda, &work[1], &
334 b[b_offset], ldb, &work[(mn << 1) + 1], info, (ftnlen)4, (ftnlen)
335 9);
336
337 /* workspace NRHS */
338
339 /* B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) */
340
341 strsm_("Left", "Upper", "No transpose", "Non-unit", rank, nrhs, &c_b36, &
342 a[a_offset], lda, &b[b_offset], ldb, (ftnlen)4, (ftnlen)5, (
343 ftnlen)12, (ftnlen)8);
344
345 i__1 = *n;
346 for (i__ = *rank + 1; i__ <= i__1; ++i__) {
347 i__2 = *nrhs;
348 for (j = 1; j <= i__2; ++j) {
349 b[i__ + j * b_dim1] = 0.f;
350 /* L30: */
351 }
352 /* L40: */
353 }
354
355 /* B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS) */
356
357 if (*rank < *n) {
358 i__1 = *rank;
359 for (i__ = 1; i__ <= i__1; ++i__) {
360 i__2 = *n - *rank + 1;
361 slatzm_("Left", &i__2, nrhs, &a[i__ + (*rank + 1) * a_dim1], lda,
362 &work[mn + i__], &b[i__ + b_dim1], &b[*rank + 1 + b_dim1],
363 ldb, &work[(mn << 1) + 1], (ftnlen)4);
364 /* L50: */
365 }
366 }
367
368 /* workspace NRHS */
369
370 /* B(1:N,1:NRHS) := P * B(1:N,1:NRHS) */
371
372 i__1 = *nrhs;
373 for (j = 1; j <= i__1; ++j) {
374 i__2 = *n;
375 for (i__ = 1; i__ <= i__2; ++i__) {
376 work[(mn << 1) + i__] = 1.f;
377 /* L60: */
378 }
379 i__2 = *n;
380 for (i__ = 1; i__ <= i__2; ++i__) {
381 if (work[(mn << 1) + i__] == 1.f) {
382 if (jpvt[i__] != i__) {
383 k = i__;
384 t1 = b[k + j * b_dim1];
385 t2 = b[jpvt[k] + j * b_dim1];
386 L70:
387 b[jpvt[k] + j * b_dim1] = t1;
388 work[(mn << 1) + k] = 0.f;
389 t1 = t2;
390 k = jpvt[k];
391 t2 = b[jpvt[k] + j * b_dim1];
392 if (jpvt[k] != i__) {
393 goto L70;
394 }
395 b[i__ + j * b_dim1] = t1;
396 work[(mn << 1) + k] = 0.f;
397 }
398 }
399 /* L80: */
400 }
401 /* L90: */
402 }
403
404 /* Undo scaling */
405
406 if (iascl == 1) {
407 slascl_("G", &c__0, &c__0, &anrm, &smlnum, n, nrhs, &b[b_offset], ldb,
408 info, (ftnlen)1);
409 slascl_("U", &c__0, &c__0, &smlnum, &anrm, rank, rank, &a[a_offset],
410 lda, info, (ftnlen)1);
411 } else if (iascl == 2) {
412 slascl_("G", &c__0, &c__0, &anrm, &bignum, n, nrhs, &b[b_offset], ldb,
413 info, (ftnlen)1);
414 slascl_("U", &c__0, &c__0, &bignum, &anrm, rank, rank, &a[a_offset],
415 lda, info, (ftnlen)1);
416 }
417 if (ibscl == 1) {
418 slascl_("G", &c__0, &c__0, &smlnum, &bnrm, n, nrhs, &b[b_offset], ldb,
419 info, (ftnlen)1);
420 } else if (ibscl == 2) {
421 slascl_("G", &c__0, &c__0, &bignum, &bnrm, n, nrhs, &b[b_offset], ldb,
422 info, (ftnlen)1);
423 }
424
425 L100:
426
427 return 0;
428
429 /* End of SGELSX */
430
431 } /* sgelsx_ */
432
433