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