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