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