1 /* linpack/cqrsl.f -- translated by f2c (version 20050501).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #ifdef __cplusplus
14 extern "C" {
15 #endif
16 #include "v3p_netlib.h"
17
18 /* Table of constant values */
19
20 static integer c__1 = 1;
21
22 /*< subroutine cqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info) >*/
cqrsl_(complex * x,integer * ldx,integer * n,integer * k,complex * qraux,complex * y,complex * qy,complex * qty,complex * b,complex * rsd,complex * xb,integer * job,integer * info)23 /* Subroutine */ int cqrsl_(complex *x, integer *ldx, integer *n, integer *k,
24 complex *qraux, complex *y, complex *qy, complex *qty, complex *b,
25 complex *rsd, complex *xb, integer *job, integer *info)
26 {
27 /* System generated locals */
28 integer x_dim1, x_offset, i__1, i__2, i__3;
29 real r__1, r__2;
30 complex q__1, q__2, q__3;
31
32 /* Builtin functions */
33 double r_imag(complex *);
34 void c_div(complex *, complex *, complex *);
35
36 /* Local variables */
37 integer i__, j;
38 complex t;
39 logical cb;
40 integer jj;
41 logical cr;
42 integer ju, kp1;
43 logical cxb, cqy;
44 complex temp;
45 logical cqty;
46 extern /* Complex */ VOID cdotc_(complex *, integer *, complex *, integer
47 *, complex *, integer *);
48 extern /* Subroutine */ int ccopy_(integer *, complex *, integer *,
49 complex *, integer *), caxpy_(integer *, complex *, complex *,
50 integer *, complex *, integer *);
51
52 /*< integer ldx,n,k,job,info >*/
53 /*< complex x(ldx,1),qraux(1),y(1),qy(1),qty(1),b(1),rsd(1),xb(1) >*/
54
55 /* cqrsl applies the output of cqrdc to compute coordinate */
56 /* transformations, projections, and least squares solutions. */
57 /* for k .le. min(n,p), let xk be the matrix */
58
59 /* xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k))) */
60
61 /* formed from columns jpvt(1), ... ,jpvt(k) of the original */
62 /* n x p matrix x that was input to cqrdc (if no pivoting was */
63 /* done, xk consists of the first k columns of x in their */
64 /* original order). cqrdc produces a factored unitary matrix q */
65 /* and an upper triangular matrix r such that */
66
67 /* xk = q * (r) */
68 /* (0) */
69
70 /* this information is contained in coded form in the arrays */
71 /* x and qraux. */
72
73 /* on entry */
74
75 /* x complex(ldx,p). */
76 /* x contains the output of cqrdc. */
77
78 /* ldx integer. */
79 /* ldx is the leading dimension of the array x. */
80
81 /* n integer. */
82 /* n is the number of rows of the matrix xk. it must */
83 /* have the same value as n in cqrdc. */
84
85 /* k integer. */
86 /* k is the number of columns of the matrix xk. k */
87 /* must nnot be greater than min(n,p), where p is the */
88 /* same as in the calling sequence to cqrdc. */
89
90 /* qraux complex(p). */
91 /* qraux contains the auxiliary output from cqrdc. */
92
93 /* y complex(n) */
94 /* y contains an n-vector that is to be manipulated */
95 /* by cqrsl. */
96
97 /* job integer. */
98 /* job specifies what is to be computed. job has */
99 /* the decimal expansion abcde, with the following */
100 /* meaning. */
101
102 /* if a.ne.0, compute qy. */
103 /* if b,c,d, or e .ne. 0, compute qty. */
104 /* if c.ne.0, compute b. */
105 /* if d.ne.0, compute rsd. */
106 /* if e.ne.0, compute xb. */
107
108 /* note that a request to compute b, rsd, or xb */
109 /* automatically triggers the computation of qty, for */
110 /* which an array must be provided in the calling */
111 /* sequence. */
112
113 /* on return */
114
115 /* qy complex(n). */
116 /* qy conntains q*y, if its computation has been */
117 /* requested. */
118
119 /* qty complex(n). */
120 /* qty contains ctrans(q)*y, if its computation has */
121 /* been requested. here ctrans(q) is the conjugate */
122 /* transpose of the matrix q. */
123
124 /* b complex(k) */
125 /* b contains the solution of the least squares problem */
126
127 /* minimize norm2(y - xk*b), */
128
129 /* if its computation has been requested. (note that */
130 /* if pivoting was requested in cqrdc, the j-th */
131 /* component of b will be associated with column jpvt(j) */
132 /* of the original matrix x that was input into cqrdc.) */
133
134 /* rsd complex(n). */
135 /* rsd contains the least squares residual y - xk*b, */
136 /* if its computation has been requested. rsd is */
137 /* also the orthogonal projection of y onto the */
138 /* orthogonal complement of the column space of xk. */
139
140 /* xb complex(n). */
141 /* xb contains the least squares approximation xk*b, */
142 /* if its computation has been requested. xb is also */
143 /* the orthogonal projection of y onto the column space */
144 /* of x. */
145
146 /* info integer. */
147 /* info is zero unless the computation of b has */
148 /* been requested and r is exactly singular. in */
149 /* this case, info is the index of the first zero */
150 /* diagonal element of r and b is left unaltered. */
151
152 /* the parameters qy, qty, b, rsd, and xb are not referenced */
153 /* if their computation is not requested and in this case */
154 /* can be replaced by dummy variables in the calling program. */
155 /* to save storage, the user may in some cases use the same */
156 /* array for different parameters in the calling sequence. a */
157 /* frequently occurring example is when one wishes to compute */
158 /* any of b, rsd, or xb and does not need y or qty. in this */
159 /* case one may identify y, qty, and one of b, rsd, or xb, while */
160 /* providing separate arrays for anything else that is to be */
161 /* computed. thus the calling sequence */
162
163 /* call cqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info) */
164
165 /* will result in the computation of b and rsd, with rsd */
166 /* overwriting y. more generally, each item in the following */
167 /* list contains groups of permissible identifications for */
168 /* a single callinng sequence. */
169
170 /* 1. (y,qty,b) (rsd) (xb) (qy) */
171
172 /* 2. (y,qty,rsd) (b) (xb) (qy) */
173
174 /* 3. (y,qty,xb) (b) (rsd) (qy) */
175
176 /* 4. (y,qy) (qty,b) (rsd) (xb) */
177
178 /* 5. (y,qy) (qty,rsd) (b) (xb) */
179
180 /* 6. (y,qy) (qty,xb) (b) (rsd) */
181
182 /* in any group the value returned in the array allocated to */
183 /* the group corresponds to the last member of the group. */
184
185 /* linpack. this version dated 08/14/78 . */
186 /* g.w. stewart, university of maryland, argonne national lab. */
187
188 /* cqrsl uses the following functions and subprograms. */
189
190 /* blas caxpy,ccopy,cdotc */
191 /* fortran abs,aimag,min0,mod,real */
192
193 /* internal variables */
194
195 /*< integer i,j,jj,ju,kp1 >*/
196 /*< complex cdotc,t,temp >*/
197 /*< logical cb,cqy,cqty,cr,cxb >*/
198
199 /*< complex zdum >*/
200 /*< real cabs1 >*/
201 /*< cabs1(zdum) = abs(real(zdum)) + abs(aimag(zdum)) >*/
202
203 /* set info flag. */
204
205 /*< info = 0 >*/
206 /* Parameter adjustments */
207 x_dim1 = *ldx;
208 x_offset = 1 + x_dim1;
209 x -= x_offset;
210 --qraux;
211 --y;
212 --qy;
213 --qty;
214 --b;
215 --rsd;
216 --xb;
217
218 /* Function Body */
219 *info = 0;
220
221 /* determine what is to be computed. */
222
223 /*< cqy = job/10000 .ne. 0 >*/
224 cqy = *job / 10000 != 0;
225 /*< cqty = mod(job,10000) .ne. 0 >*/
226 cqty = *job % 10000 != 0;
227 /*< cb = mod(job,1000)/100 .ne. 0 >*/
228 cb = *job % 1000 / 100 != 0;
229 /*< cr = mod(job,100)/10 .ne. 0 >*/
230 cr = *job % 100 / 10 != 0;
231 /*< cxb = mod(job,10) .ne. 0 >*/
232 cxb = *job % 10 != 0;
233 /*< ju = min0(k,n-1) >*/
234 /* Computing MIN */
235 i__1 = *k, i__2 = *n - 1;
236 ju = min(i__1,i__2);
237
238 /* special action when n=1. */
239
240 /*< if (ju .ne. 0) go to 40 >*/
241 if (ju != 0) {
242 goto L40;
243 }
244 /*< if (cqy) qy(1) = y(1) >*/
245 if (cqy) {
246 qy[1].r = y[1].r, qy[1].i = y[1].i;
247 }
248 /*< if (cqty) qty(1) = y(1) >*/
249 if (cqty) {
250 qty[1].r = y[1].r, qty[1].i = y[1].i;
251 }
252 /*< if (cxb) xb(1) = y(1) >*/
253 if (cxb) {
254 xb[1].r = y[1].r, xb[1].i = y[1].i;
255 }
256 /*< if (.not.cb) go to 30 >*/
257 if (! cb) {
258 goto L30;
259 }
260 /*< if (cabs1(x(1,1)) .ne. 0.0e0) go to 10 >*/
261 i__1 = x_dim1 + 1;
262 if ((r__1 = x[i__1].r, dabs(r__1)) + (r__2 = r_imag(&x[x_dim1 + 1]), dabs(
263 r__2)) != (float)0.) {
264 goto L10;
265 }
266 /*< info = 1 >*/
267 *info = 1;
268 /*< go to 20 >*/
269 goto L20;
270 /*< 10 continue >*/
271 L10:
272 /*< b(1) = y(1)/x(1,1) >*/
273 c_div(&q__1, &y[1], &x[x_dim1 + 1]);
274 b[1].r = q__1.r, b[1].i = q__1.i;
275 /*< 20 continue >*/
276 L20:
277 /*< 30 continue >*/
278 L30:
279 /*< if (cr) rsd(1) = (0.0e0,0.0e0) >*/
280 if (cr) {
281 rsd[1].r = (float)0., rsd[1].i = (float)0.;
282 }
283 /*< go to 250 >*/
284 goto L250;
285 /*< 40 continue >*/
286 L40:
287
288 /* set up to compute qy or qty. */
289
290 /*< if (cqy) call ccopy(n,y,1,qy,1) >*/
291 if (cqy) {
292 ccopy_(n, &y[1], &c__1, &qy[1], &c__1);
293 }
294 /*< if (cqty) call ccopy(n,y,1,qty,1) >*/
295 if (cqty) {
296 ccopy_(n, &y[1], &c__1, &qty[1], &c__1);
297 }
298 /*< if (.not.cqy) go to 70 >*/
299 if (! cqy) {
300 goto L70;
301 }
302
303 /* compute qy. */
304
305 /*< do 60 jj = 1, ju >*/
306 i__1 = ju;
307 for (jj = 1; jj <= i__1; ++jj) {
308 /*< j = ju - jj + 1 >*/
309 j = ju - jj + 1;
310 /*< if (cabs1(qraux(j)) .eq. 0.0e0) go to 50 >*/
311 i__2 = j;
312 if ((r__1 = qraux[i__2].r, dabs(r__1)) + (r__2 = r_imag(&qraux[j]),
313 dabs(r__2)) == (float)0.) {
314 goto L50;
315 }
316 /*< temp = x(j,j) >*/
317 i__2 = j + j * x_dim1;
318 temp.r = x[i__2].r, temp.i = x[i__2].i;
319 /*< x(j,j) = qraux(j) >*/
320 i__2 = j + j * x_dim1;
321 i__3 = j;
322 x[i__2].r = qraux[i__3].r, x[i__2].i = qraux[i__3].i;
323 /*< t = -cdotc(n-j+1,x(j,j),1,qy(j),1)/x(j,j) >*/
324 i__2 = *n - j + 1;
325 cdotc_(&q__3, &i__2, &x[j + j * x_dim1], &c__1, &qy[j], &c__1);
326 q__2.r = -q__3.r, q__2.i = -q__3.i;
327 c_div(&q__1, &q__2, &x[j + j * x_dim1]);
328 t.r = q__1.r, t.i = q__1.i;
329 /*< call caxpy(n-j+1,t,x(j,j),1,qy(j),1) >*/
330 i__2 = *n - j + 1;
331 caxpy_(&i__2, &t, &x[j + j * x_dim1], &c__1, &qy[j], &c__1);
332 /*< x(j,j) = temp >*/
333 i__2 = j + j * x_dim1;
334 x[i__2].r = temp.r, x[i__2].i = temp.i;
335 /*< 50 continue >*/
336 L50:
337 /*< 60 continue >*/
338 /* L60: */
339 ;
340 }
341 /*< 70 continue >*/
342 L70:
343 /*< if (.not.cqty) go to 100 >*/
344 if (! cqty) {
345 goto L100;
346 }
347
348 /* compute ctrans(q)*y. */
349
350 /*< do 90 j = 1, ju >*/
351 i__1 = ju;
352 for (j = 1; j <= i__1; ++j) {
353 /*< if (cabs1(qraux(j)) .eq. 0.0e0) go to 80 >*/
354 i__2 = j;
355 if ((r__1 = qraux[i__2].r, dabs(r__1)) + (r__2 = r_imag(&qraux[j]),
356 dabs(r__2)) == (float)0.) {
357 goto L80;
358 }
359 /*< temp = x(j,j) >*/
360 i__2 = j + j * x_dim1;
361 temp.r = x[i__2].r, temp.i = x[i__2].i;
362 /*< x(j,j) = qraux(j) >*/
363 i__2 = j + j * x_dim1;
364 i__3 = j;
365 x[i__2].r = qraux[i__3].r, x[i__2].i = qraux[i__3].i;
366 /*< t = -cdotc(n-j+1,x(j,j),1,qty(j),1)/x(j,j) >*/
367 i__2 = *n - j + 1;
368 cdotc_(&q__3, &i__2, &x[j + j * x_dim1], &c__1, &qty[j], &c__1);
369 q__2.r = -q__3.r, q__2.i = -q__3.i;
370 c_div(&q__1, &q__2, &x[j + j * x_dim1]);
371 t.r = q__1.r, t.i = q__1.i;
372 /*< call caxpy(n-j+1,t,x(j,j),1,qty(j),1) >*/
373 i__2 = *n - j + 1;
374 caxpy_(&i__2, &t, &x[j + j * x_dim1], &c__1, &qty[j], &c__1);
375 /*< x(j,j) = temp >*/
376 i__2 = j + j * x_dim1;
377 x[i__2].r = temp.r, x[i__2].i = temp.i;
378 /*< 80 continue >*/
379 L80:
380 /*< 90 continue >*/
381 /* L90: */
382 ;
383 }
384 /*< 100 continue >*/
385 L100:
386
387 /* set up to compute b, rsd, or xb. */
388
389 /*< if (cb) call ccopy(k,qty,1,b,1) >*/
390 if (cb) {
391 ccopy_(k, &qty[1], &c__1, &b[1], &c__1);
392 }
393 /*< kp1 = k + 1 >*/
394 kp1 = *k + 1;
395 /*< if (cxb) call ccopy(k,qty,1,xb,1) >*/
396 if (cxb) {
397 ccopy_(k, &qty[1], &c__1, &xb[1], &c__1);
398 }
399 /*< if (cr .and. k .lt. n) call ccopy(n-k,qty(kp1),1,rsd(kp1),1) >*/
400 if (cr && *k < *n) {
401 i__1 = *n - *k;
402 ccopy_(&i__1, &qty[kp1], &c__1, &rsd[kp1], &c__1);
403 }
404 /*< if (.not.cxb .or. kp1 .gt. n) go to 120 >*/
405 if (! cxb || kp1 > *n) {
406 goto L120;
407 }
408 /*< do 110 i = kp1, n >*/
409 i__1 = *n;
410 for (i__ = kp1; i__ <= i__1; ++i__) {
411 /*< xb(i) = (0.0e0,0.0e0) >*/
412 i__2 = i__;
413 xb[i__2].r = (float)0., xb[i__2].i = (float)0.;
414 /*< 110 continue >*/
415 /* L110: */
416 }
417 /*< 120 continue >*/
418 L120:
419 /*< if (.not.cr) go to 140 >*/
420 if (! cr) {
421 goto L140;
422 }
423 /*< do 130 i = 1, k >*/
424 i__1 = *k;
425 for (i__ = 1; i__ <= i__1; ++i__) {
426 /*< rsd(i) = (0.0e0,0.0e0) >*/
427 i__2 = i__;
428 rsd[i__2].r = (float)0., rsd[i__2].i = (float)0.;
429 /*< 130 continue >*/
430 /* L130: */
431 }
432 /*< 140 continue >*/
433 L140:
434 /*< if (.not.cb) go to 190 >*/
435 if (! cb) {
436 goto L190;
437 }
438
439 /* compute b. */
440
441 /*< do 170 jj = 1, k >*/
442 i__1 = *k;
443 for (jj = 1; jj <= i__1; ++jj) {
444 /*< j = k - jj + 1 >*/
445 j = *k - jj + 1;
446 /*< if (cabs1(x(j,j)) .ne. 0.0e0) go to 150 >*/
447 i__2 = j + j * x_dim1;
448 if ((r__1 = x[i__2].r, dabs(r__1)) + (r__2 = r_imag(&x[j + j * x_dim1]
449 ), dabs(r__2)) != (float)0.) {
450 goto L150;
451 }
452 /*< info = j >*/
453 *info = j;
454 /* ......exit */
455 /*< go to 180 >*/
456 goto L180;
457 /*< 150 continue >*/
458 L150:
459 /*< b(j) = b(j)/x(j,j) >*/
460 i__2 = j;
461 c_div(&q__1, &b[j], &x[j + j * x_dim1]);
462 b[i__2].r = q__1.r, b[i__2].i = q__1.i;
463 /*< if (j .eq. 1) go to 160 >*/
464 if (j == 1) {
465 goto L160;
466 }
467 /*< t = -b(j) >*/
468 i__2 = j;
469 q__1.r = -b[i__2].r, q__1.i = -b[i__2].i;
470 t.r = q__1.r, t.i = q__1.i;
471 /*< call caxpy(j-1,t,x(1,j),1,b,1) >*/
472 i__2 = j - 1;
473 caxpy_(&i__2, &t, &x[j * x_dim1 + 1], &c__1, &b[1], &c__1);
474 /*< 160 continue >*/
475 L160:
476 /*< 170 continue >*/
477 /* L170: */
478 ;
479 }
480 /*< 180 continue >*/
481 L180:
482 /*< 190 continue >*/
483 L190:
484 /*< if (.not.cr .and. .not.cxb) go to 240 >*/
485 if (! cr && ! cxb) {
486 goto L240;
487 }
488
489 /* compute rsd or xb as required. */
490
491 /*< do 230 jj = 1, ju >*/
492 i__1 = ju;
493 for (jj = 1; jj <= i__1; ++jj) {
494 /*< j = ju - jj + 1 >*/
495 j = ju - jj + 1;
496 /*< if (cabs1(qraux(j)) .eq. 0.0e0) go to 220 >*/
497 i__2 = j;
498 if ((r__1 = qraux[i__2].r, dabs(r__1)) + (r__2 = r_imag(&qraux[j]),
499 dabs(r__2)) == (float)0.) {
500 goto L220;
501 }
502 /*< temp = x(j,j) >*/
503 i__2 = j + j * x_dim1;
504 temp.r = x[i__2].r, temp.i = x[i__2].i;
505 /*< x(j,j) = qraux(j) >*/
506 i__2 = j + j * x_dim1;
507 i__3 = j;
508 x[i__2].r = qraux[i__3].r, x[i__2].i = qraux[i__3].i;
509 /*< if (.not.cr) go to 200 >*/
510 if (! cr) {
511 goto L200;
512 }
513 /*< t = -cdotc(n-j+1,x(j,j),1,rsd(j),1)/x(j,j) >*/
514 i__2 = *n - j + 1;
515 cdotc_(&q__3, &i__2, &x[j + j * x_dim1], &c__1, &rsd[j], &c__1);
516 q__2.r = -q__3.r, q__2.i = -q__3.i;
517 c_div(&q__1, &q__2, &x[j + j * x_dim1]);
518 t.r = q__1.r, t.i = q__1.i;
519 /*< call caxpy(n-j+1,t,x(j,j),1,rsd(j),1) >*/
520 i__2 = *n - j + 1;
521 caxpy_(&i__2, &t, &x[j + j * x_dim1], &c__1, &rsd[j], &c__1);
522 /*< 200 continue >*/
523 L200:
524 /*< if (.not.cxb) go to 210 >*/
525 if (! cxb) {
526 goto L210;
527 }
528 /*< t = -cdotc(n-j+1,x(j,j),1,xb(j),1)/x(j,j) >*/
529 i__2 = *n - j + 1;
530 cdotc_(&q__3, &i__2, &x[j + j * x_dim1], &c__1, &xb[j], &c__1);
531 q__2.r = -q__3.r, q__2.i = -q__3.i;
532 c_div(&q__1, &q__2, &x[j + j * x_dim1]);
533 t.r = q__1.r, t.i = q__1.i;
534 /*< call caxpy(n-j+1,t,x(j,j),1,xb(j),1) >*/
535 i__2 = *n - j + 1;
536 caxpy_(&i__2, &t, &x[j + j * x_dim1], &c__1, &xb[j], &c__1);
537 /*< 210 continue >*/
538 L210:
539 /*< x(j,j) = temp >*/
540 i__2 = j + j * x_dim1;
541 x[i__2].r = temp.r, x[i__2].i = temp.i;
542 /*< 220 continue >*/
543 L220:
544 /*< 230 continue >*/
545 /* L230: */
546 ;
547 }
548 /*< 240 continue >*/
549 L240:
550 /*< 250 continue >*/
551 L250:
552 /*< return >*/
553 return 0;
554 /*< end >*/
555 } /* cqrsl_ */
556
557 #ifdef __cplusplus
558 }
559 #endif
560