1 /* linpack/zqrdc.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 static doublecomplex c_b28 = {1.,0.};
22
23 /*< subroutine zqrdc(x,ldx,n,p,qraux,jpvt,work,job) >*/
zqrdc_(doublecomplex * x,integer * ldx,integer * n,integer * p,doublecomplex * qraux,integer * jpvt,doublecomplex * work,integer * job)24 /* Subroutine */ int zqrdc_(doublecomplex *x, integer *ldx, integer *n,
25 integer *p, doublecomplex *qraux, integer *jpvt, doublecomplex *work,
26 integer *job)
27 {
28 /* System generated locals */
29 integer x_dim1, x_offset, i__1, i__2, i__3, i__4;
30 doublereal d__1, d__2, d__3, d__4;
31 doublecomplex z__1, z__2, z__3;
32
33 /* Builtin functions */
34 double z_abs(doublecomplex *);
35 void z_div(doublecomplex *, doublecomplex *, doublecomplex *), z_sqrt(
36 doublecomplex *, doublecomplex *);
37
38 /* Local variables */
39 integer j, l;
40 doublecomplex t;
41 integer jj, jp, pl, pu;
42 doublereal tt;
43 integer lp1, lup;
44 logical negj;
45 integer maxj;
46 extern /* Subroutine */ int zscal_(integer *, doublecomplex *,
47 doublecomplex *, integer *);
48 extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *,
49 doublecomplex *, integer *, doublecomplex *, integer *);
50 logical swapj;
51 doublecomplex nrmxl;
52 extern /* Subroutine */ int zswap_(integer *, doublecomplex *, integer *,
53 doublecomplex *, integer *), zaxpy_(integer *, doublecomplex *,
54 doublecomplex *, integer *, doublecomplex *, integer *);
55 extern doublereal dznrm2_(integer *, doublecomplex *, integer *);
56 doublereal maxnrm;
57
58 /*< integer ldx,n,p,job >*/
59 /*< integer jpvt(1) >*/
60 /*< complex*16 x(ldx,1),qraux(1),work(1) >*/
61
62 /* zqrdc uses householder transformations to compute the qr */
63 /* factorization of an n by p matrix x. column pivoting */
64 /* based on the 2-norms of the reduced columns may be */
65 /* performed at the users option. */
66
67 /* on entry */
68
69 /* x complex*16(ldx,p), where ldx .ge. n. */
70 /* x contains the matrix whose decomposition is to be */
71 /* computed. */
72
73 /* ldx integer. */
74 /* ldx is the leading dimension of the array x. */
75
76 /* n integer. */
77 /* n is the number of rows of the matrix x. */
78
79 /* p integer. */
80 /* p is the number of columns of the matrix x. */
81
82 /* jpvt integer(p). */
83 /* jpvt contains integers that control the selection */
84 /* of the pivot columns. the k-th column x(k) of x */
85 /* is placed in one of three classes according to the */
86 /* value of jpvt(k). */
87
88 /* if jpvt(k) .gt. 0, then x(k) is an initial */
89 /* column. */
90
91 /* if jpvt(k) .eq. 0, then x(k) is a free column. */
92
93 /* if jpvt(k) .lt. 0, then x(k) is a final column. */
94
95 /* before the decomposition is computed, initial columns */
96 /* are moved to the beginning of the array x and final */
97 /* columns to the end. both initial and final columns */
98 /* are frozen in place during the computation and only */
99 /* free columns are moved. at the k-th stage of the */
100 /* reduction, if x(k) is occupied by a free column */
101 /* it is interchanged with the free column of largest */
102 /* reduced norm. jpvt is not referenced if */
103 /* job .eq. 0. */
104
105 /* work complex*16(p). */
106 /* work is a work array. work is not referenced if */
107 /* job .eq. 0. */
108
109 /* job integer. */
110 /* job is an integer that initiates column pivoting. */
111 /* if job .eq. 0, no pivoting is done. */
112 /* if job .ne. 0, pivoting is done. */
113
114 /* on return */
115
116 /* x x contains in its upper triangle the upper */
117 /* triangular matrix r of the qr factorization. */
118 /* below its diagonal x contains information from */
119 /* which the unitary part of the decomposition */
120 /* can be recovered. note that if pivoting has */
121 /* been requested, the decomposition is not that */
122 /* of the original matrix x but that of x */
123 /* with its columns permuted as described by jpvt. */
124
125 /* qraux complex*16(p). */
126 /* qraux contains further information required to recover */
127 /* the unitary part of the decomposition. */
128
129 /* jpvt jpvt(k) contains the index of the column of the */
130 /* original matrix that has been interchanged into */
131 /* the k-th column, if pivoting was requested. */
132
133 /* linpack. this version dated 08/14/78 . */
134 /* g.w. stewart, university of maryland, argonne national lab. */
135
136 /* zqrdc uses the following functions and subprograms. */
137
138 /* blas zaxpy,zdotc,zscal,zswap,dznrm2 */
139 /* fortran dabs,dmax1,cdabs,dcmplx,cdsqrt,min0 */
140
141 /* internal variables */
142
143 /*< integer j,jp,l,lp1,lup,maxj,pl,pu >*/
144 /*< double precision maxnrm,dznrm2,tt >*/
145 /*< complex*16 zdotc,nrmxl,t >*/
146 /*< logical negj,swapj >*/
147
148 /*< complex*16 csign,zdum,zdum1,zdum2 >*/
149 /*< double precision cabs1 >*/
150 /*< double precision dreal,dimag >*/
151 /*< complex*16 zdumr,zdumi >*/
152 /*< dreal(zdumr) = zdumr >*/
153 /*< dimag(zdumi) = (0.0d0,-1.0d0)*zdumi >*/
154 /*< csign(zdum1,zdum2) = cdabs(zdum1)*(zdum2/cdabs(zdum2)) >*/
155 /*< cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum)) >*/
156
157 /*< pl = 1 >*/
158 /* Parameter adjustments */
159 x_dim1 = *ldx;
160 x_offset = 1 + x_dim1;
161 x -= x_offset;
162 --qraux;
163 --jpvt;
164 --work;
165
166 /* Function Body */
167 pl = 1;
168 /*< pu = 0 >*/
169 pu = 0;
170 /*< if (job .eq. 0) go to 60 >*/
171 if (*job == 0) {
172 goto L60;
173 }
174
175 /* pivoting has been requested. rearrange the columns */
176 /* according to jpvt. */
177
178 /*< do 20 j = 1, p >*/
179 i__1 = *p;
180 for (j = 1; j <= i__1; ++j) {
181 /*< swapj = jpvt(j) .gt. 0 >*/
182 swapj = jpvt[j] > 0;
183 /*< negj = jpvt(j) .lt. 0 >*/
184 negj = jpvt[j] < 0;
185 /*< jpvt(j) = j >*/
186 jpvt[j] = j;
187 /*< if (negj) jpvt(j) = -j >*/
188 if (negj) {
189 jpvt[j] = -j;
190 }
191 /*< if (.not.swapj) go to 10 >*/
192 if (! swapj) {
193 goto L10;
194 }
195 /*< if (j .ne. pl) call zswap(n,x(1,pl),1,x(1,j),1) >*/
196 if (j != pl) {
197 zswap_(n, &x[pl * x_dim1 + 1], &c__1, &x[j * x_dim1 + 1], &c__1);
198 }
199 /*< jpvt(j) = jpvt(pl) >*/
200 jpvt[j] = jpvt[pl];
201 /*< jpvt(pl) = j >*/
202 jpvt[pl] = j;
203 /*< pl = pl + 1 >*/
204 ++pl;
205 /*< 10 continue >*/
206 L10:
207 /*< 20 continue >*/
208 /* L20: */
209 ;
210 }
211 /*< pu = p >*/
212 pu = *p;
213 /*< do 50 jj = 1, p >*/
214 i__1 = *p;
215 for (jj = 1; jj <= i__1; ++jj) {
216 /*< j = p - jj + 1 >*/
217 j = *p - jj + 1;
218 /*< if (jpvt(j) .ge. 0) go to 40 >*/
219 if (jpvt[j] >= 0) {
220 goto L40;
221 }
222 /*< jpvt(j) = -jpvt(j) >*/
223 jpvt[j] = -jpvt[j];
224 /*< if (j .eq. pu) go to 30 >*/
225 if (j == pu) {
226 goto L30;
227 }
228 /*< call zswap(n,x(1,pu),1,x(1,j),1) >*/
229 zswap_(n, &x[pu * x_dim1 + 1], &c__1, &x[j * x_dim1 + 1], &c__1);
230 /*< jp = jpvt(pu) >*/
231 jp = jpvt[pu];
232 /*< jpvt(pu) = jpvt(j) >*/
233 jpvt[pu] = jpvt[j];
234 /*< jpvt(j) = jp >*/
235 jpvt[j] = jp;
236 /*< 30 continue >*/
237 L30:
238 /*< pu = pu - 1 >*/
239 --pu;
240 /*< 40 continue >*/
241 L40:
242 /*< 50 continue >*/
243 /* L50: */
244 ;
245 }
246 /*< 60 continue >*/
247 L60:
248
249 /* compute the norms of the free columns. */
250
251 /*< if (pu .lt. pl) go to 80 >*/
252 if (pu < pl) {
253 goto L80;
254 }
255 /*< do 70 j = pl, pu >*/
256 i__1 = pu;
257 for (j = pl; j <= i__1; ++j) {
258 /*< qraux(j) = dcmplx(dznrm2(n,x(1,j),1),0.0d0) >*/
259 i__2 = j;
260 d__1 = dznrm2_(n, &x[j * x_dim1 + 1], &c__1);
261 z__1.r = d__1, z__1.i = 0.;
262 qraux[i__2].r = z__1.r, qraux[i__2].i = z__1.i;
263 /*< work(j) = qraux(j) >*/
264 i__2 = j;
265 i__3 = j;
266 work[i__2].r = qraux[i__3].r, work[i__2].i = qraux[i__3].i;
267 /*< 70 continue >*/
268 /* L70: */
269 }
270 /*< 80 continue >*/
271 L80:
272
273 /* perform the householder reduction of x. */
274
275 /*< lup = min0(n,p) >*/
276 lup = min(*n,*p);
277 /*< do 200 l = 1, lup >*/
278 i__1 = lup;
279 for (l = 1; l <= i__1; ++l) {
280 /*< if (l .lt. pl .or. l .ge. pu) go to 120 >*/
281 if (l < pl || l >= pu) {
282 goto L120;
283 }
284
285 /* locate the column of largest norm and bring it */
286 /* into the pivot position. */
287
288 /*< maxnrm = 0.0d0 >*/
289 maxnrm = 0.;
290 /*< maxj = l >*/
291 maxj = l;
292 /*< do 100 j = l, pu >*/
293 i__2 = pu;
294 for (j = l; j <= i__2; ++j) {
295 /*< if (dreal(qraux(j)) .le. maxnrm) go to 90 >*/
296 i__3 = j;
297 if (qraux[i__3].r <= maxnrm) {
298 goto L90;
299 }
300 /*< maxnrm = dreal(qraux(j)) >*/
301 i__3 = j;
302 maxnrm = qraux[i__3].r;
303 /*< maxj = j >*/
304 maxj = j;
305 /*< 90 continue >*/
306 L90:
307 /*< 100 continue >*/
308 /* L100: */
309 ;
310 }
311 /*< if (maxj .eq. l) go to 110 >*/
312 if (maxj == l) {
313 goto L110;
314 }
315 /*< call zswap(n,x(1,l),1,x(1,maxj),1) >*/
316 zswap_(n, &x[l * x_dim1 + 1], &c__1, &x[maxj * x_dim1 + 1], &c__1);
317 /*< qraux(maxj) = qraux(l) >*/
318 i__2 = maxj;
319 i__3 = l;
320 qraux[i__2].r = qraux[i__3].r, qraux[i__2].i = qraux[i__3].i;
321 /*< work(maxj) = work(l) >*/
322 i__2 = maxj;
323 i__3 = l;
324 work[i__2].r = work[i__3].r, work[i__2].i = work[i__3].i;
325 /*< jp = jpvt(maxj) >*/
326 jp = jpvt[maxj];
327 /*< jpvt(maxj) = jpvt(l) >*/
328 jpvt[maxj] = jpvt[l];
329 /*< jpvt(l) = jp >*/
330 jpvt[l] = jp;
331 /*< 110 continue >*/
332 L110:
333 /*< 120 continue >*/
334 L120:
335 /*< qraux(l) = (0.0d0,0.0d0) >*/
336 i__2 = l;
337 qraux[i__2].r = 0., qraux[i__2].i = 0.;
338 /*< if (l .eq. n) go to 190 >*/
339 if (l == *n) {
340 goto L190;
341 }
342
343 /* compute the householder transformation for column l. */
344
345 /*< nrmxl = dcmplx(dznrm2(n-l+1,x(l,l),1),0.0d0) >*/
346 i__2 = *n - l + 1;
347 d__1 = dznrm2_(&i__2, &x[l + l * x_dim1], &c__1);
348 z__1.r = d__1, z__1.i = 0.;
349 nrmxl.r = z__1.r, nrmxl.i = z__1.i;
350 /*< if (cabs1(nrmxl) .eq. 0.0d0) go to 180 >*/
351 z__1.r = nrmxl.r * 0. - nrmxl.i * -1., z__1.i = nrmxl.i * 0. +
352 nrmxl.r * -1.;
353 if ((d__1 = nrmxl.r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) == 0.) {
354 goto L180;
355 }
356 /*< >*/
357 i__2 = l + l * x_dim1;
358 i__3 = l + l * x_dim1;
359 z__1.r = x[i__3].r * 0. - x[i__3].i * -1., z__1.i = x[i__3].i * 0. +
360 x[i__3].r * -1.;
361 if ((d__1 = x[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) != 0.)
362 {
363 d__3 = z_abs(&nrmxl);
364 i__4 = l + l * x_dim1;
365 d__4 = z_abs(&x[l + l * x_dim1]);
366 z__3.r = x[i__4].r / d__4, z__3.i = x[i__4].i / d__4;
367 z__2.r = d__3 * z__3.r, z__2.i = d__3 * z__3.i;
368 nrmxl.r = z__2.r, nrmxl.i = z__2.i;
369 }
370 /*< call zscal(n-l+1,(1.0d0,0.0d0)/nrmxl,x(l,l),1) >*/
371 i__2 = *n - l + 1;
372 z_div(&z__1, &c_b28, &nrmxl);
373 zscal_(&i__2, &z__1, &x[l + l * x_dim1], &c__1);
374 /*< x(l,l) = (1.0d0,0.0d0) + x(l,l) >*/
375 i__2 = l + l * x_dim1;
376 i__3 = l + l * x_dim1;
377 z__1.r = x[i__3].r + 1., z__1.i = x[i__3].i + 0.;
378 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
379
380 /* apply the transformation to the remaining columns, */
381 /* updating the norms. */
382
383 /*< lp1 = l + 1 >*/
384 lp1 = l + 1;
385 /*< if (p .lt. lp1) go to 170 >*/
386 if (*p < lp1) {
387 goto L170;
388 }
389 /*< do 160 j = lp1, p >*/
390 i__2 = *p;
391 for (j = lp1; j <= i__2; ++j) {
392 /*< t = -zdotc(n-l+1,x(l,l),1,x(l,j),1)/x(l,l) >*/
393 i__3 = *n - l + 1;
394 zdotc_(&z__3, &i__3, &x[l + l * x_dim1], &c__1, &x[l + j * x_dim1]
395 , &c__1);
396 z__2.r = -z__3.r, z__2.i = -z__3.i;
397 z_div(&z__1, &z__2, &x[l + l * x_dim1]);
398 t.r = z__1.r, t.i = z__1.i;
399 /*< call zaxpy(n-l+1,t,x(l,l),1,x(l,j),1) >*/
400 i__3 = *n - l + 1;
401 zaxpy_(&i__3, &t, &x[l + l * x_dim1], &c__1, &x[l + j * x_dim1], &
402 c__1);
403 /*< if (j .lt. pl .or. j .gt. pu) go to 150 >*/
404 if (j < pl || j > pu) {
405 goto L150;
406 }
407 /*< if (cabs1(qraux(j)) .eq. 0.0d0) go to 150 >*/
408 i__3 = j;
409 i__4 = j;
410 z__1.r = qraux[i__4].r * 0. - qraux[i__4].i * -1., z__1.i = qraux[
411 i__4].i * 0. + qraux[i__4].r * -1.;
412 if ((d__1 = qraux[i__3].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2))
413 == 0.) {
414 goto L150;
415 }
416 /*< tt = 1.0d0 - (cdabs(x(l,j))/dreal(qraux(j)))**2 >*/
417 i__3 = j;
418 /* Computing 2nd power */
419 d__1 = z_abs(&x[l + j * x_dim1]) / qraux[i__3].r;
420 tt = 1. - d__1 * d__1;
421 /*< tt = dmax1(tt,0.0d0) >*/
422 tt = max(tt,0.);
423 /*< t = dcmplx(tt,0.0d0) >*/
424 z__1.r = tt, z__1.i = 0.;
425 t.r = z__1.r, t.i = z__1.i;
426 /*< >*/
427 i__3 = j;
428 i__4 = j;
429 /* Computing 2nd power */
430 d__1 = qraux[i__3].r / work[i__4].r;
431 tt = tt * .05 * (d__1 * d__1) + 1.;
432 /*< if (tt .eq. 1.0d0) go to 130 >*/
433 if (tt == 1.) {
434 goto L130;
435 }
436 /*< qraux(j) = qraux(j)*cdsqrt(t) >*/
437 i__3 = j;
438 i__4 = j;
439 z_sqrt(&z__2, &t);
440 z__1.r = qraux[i__4].r * z__2.r - qraux[i__4].i * z__2.i, z__1.i =
441 qraux[i__4].r * z__2.i + qraux[i__4].i * z__2.r;
442 qraux[i__3].r = z__1.r, qraux[i__3].i = z__1.i;
443 /*< go to 140 >*/
444 goto L140;
445 /*< 130 continue >*/
446 L130:
447 /*< qraux(j) = dcmplx(dznrm2(n-l,x(l+1,j),1),0.0d0) >*/
448 i__3 = j;
449 i__4 = *n - l;
450 d__1 = dznrm2_(&i__4, &x[l + 1 + j * x_dim1], &c__1);
451 z__1.r = d__1, z__1.i = 0.;
452 qraux[i__3].r = z__1.r, qraux[i__3].i = z__1.i;
453 /*< work(j) = qraux(j) >*/
454 i__3 = j;
455 i__4 = j;
456 work[i__3].r = qraux[i__4].r, work[i__3].i = qraux[i__4].i;
457 /*< 140 continue >*/
458 L140:
459 /*< 150 continue >*/
460 L150:
461 /*< 160 continue >*/
462 /* L160: */
463 ;
464 }
465 /*< 170 continue >*/
466 L170:
467
468 /* save the transformation. */
469
470 /*< qraux(l) = x(l,l) >*/
471 i__2 = l;
472 i__3 = l + l * x_dim1;
473 qraux[i__2].r = x[i__3].r, qraux[i__2].i = x[i__3].i;
474 /*< x(l,l) = -nrmxl >*/
475 i__2 = l + l * x_dim1;
476 z__1.r = -nrmxl.r, z__1.i = -nrmxl.i;
477 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
478 /*< 180 continue >*/
479 L180:
480 /*< 190 continue >*/
481 L190:
482 /*< 200 continue >*/
483 /* L200: */
484 ;
485 }
486 /*< return >*/
487 return 0;
488 /*< end >*/
489 } /* zqrdc_ */
490
491 #ifdef __cplusplus
492 }
493 #endif
494