source: orange/source/orange/linpack/dqrsl.c @ 8069:5e19a5dcaac0

Revision 8069:5e19a5dcaac0, 10.6 KB checked in by ales_erjavec <ales.erjavec@…>, 3 years ago (diff)
Line 
1/* dqrsl.f -- translated by f2c (version 20090411).
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#include "linpack.h"
14//#include "f2c.h"
15
16/* Table of constant values */
17
18static integer c__1 = 1;
19
20/* Subroutine */ int dqrsl_(doublereal *x, integer *ldx, integer *n, integer *
21    k, doublereal *qraux, doublereal *y, doublereal *qy, doublereal *qty,
22    doublereal *b, doublereal *rsd, doublereal *xb, integer *job, integer
23    *info)
24{
25    /* System generated locals */
26    integer x_dim1, x_offset, i__1, i__2;
27
28    /* Local variables */
29    static integer i__, j;
30    static doublereal t;
31    static logical cb;
32    static integer jj;
33    static logical cr;
34    static integer ju, kp1;
35    static logical cxb, cqy;
36    extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
37        integer *);
38    static doublereal temp;
39    static logical cqty;
40    extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
41        doublereal *, integer *), daxpy_(integer *, doublereal *,
42        doublereal *, integer *, doublereal *, integer *);
43
44
45/*     dqrsl applies the output of dqrdc to compute coordinate */
46/*     transformations, projections, and least squares solutions. */
47/*     for k .le. min(n,p), let xk be the matrix */
48
49/*            xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k))) */
50
51/*     formed from columnns jpvt(1), ... ,jpvt(k) of the original */
52/*     n x p matrix x that was input to dqrdc (if no pivoting was */
53/*     done, xk consists of the first k columns of x in their */
54/*     original order).  dqrdc produces a factored orthogonal matrix q */
55/*     and an upper triangular matrix r such that */
56
57/*              xk = q * (r) */
58/*                       (0) */
59
60/*     this information is contained in coded form in the arrays */
61/*     x and qraux. */
62
63/*     on entry */
64
65/*        x      double precision(ldx,p). */
66/*               x contains the output of dqrdc. */
67
68/*        ldx    integer. */
69/*               ldx is the leading dimension of the array x. */
70
71/*        n      integer. */
72/*               n is the number of rows of the matrix xk.  it must */
73/*               have the same value as n in dqrdc. */
74
75/*        k      integer. */
76/*               k is the number of columns of the matrix xk.  k */
77/*               must nnot be greater than min(n,p), where p is the */
78/*               same as in the calling sequence to dqrdc. */
79
80/*        qraux  double precision(p). */
81/*               qraux contains the auxiliary output from dqrdc. */
82
83/*        y      double precision(n) */
84/*               y contains an n-vector that is to be manipulated */
85/*               by dqrsl. */
86
87/*        job    integer. */
88/*               job specifies what is to be computed.  job has */
89/*               the decimal expansion abcde, with the following */
90/*               meaning. */
91
92/*                    if a.ne.0, compute qy. */
93/*                    if b,c,d, or e .ne. 0, compute qty. */
94/*                    if c.ne.0, compute b. */
95/*                    if d.ne.0, compute rsd. */
96/*                    if e.ne.0, compute xb. */
97
98/*               note that a request to compute b, rsd, or xb */
99/*               automatically triggers the computation of qty, for */
100/*               which an array must be provided in the calling */
101/*               sequence. */
102
103/*     on return */
104
105/*        qy     double precision(n). */
106/*               qy conntains q*y, if its computation has been */
107/*               requested. */
108
109/*        qty    double precision(n). */
110/*               qty contains trans(q)*y, if its computation has */
111/*               been requested.  here trans(q) is the */
112/*               transpose of the matrix q. */
113
114/*        b      double precision(k) */
115/*               b contains the solution of the least squares problem */
116
117/*                    minimize norm2(y - xk*b), */
118
119/*               if its computation has been requested.  (note that */
120/*               if pivoting was requested in dqrdc, the j-th */
121/*               component of b will be associated with column jpvt(j) */
122/*               of the original matrix x that was input into dqrdc.) */
123
124/*        rsd    double precision(n). */
125/*               rsd contains the least squares residual y - xk*b, */
126/*               if its computation has been requested.  rsd is */
127/*               also the orthogonal projection of y onto the */
128/*               orthogonal complement of the column space of xk. */
129
130/*        xb     double precision(n). */
131/*               xb contains the least squares approximation xk*b, */
132/*               if its computation has been requested.  xb is also */
133/*               the orthogonal projection of y onto the column space */
134/*               of x. */
135
136/*        info   integer. */
137/*               info is zero unless the computation of b has */
138/*               been requested and r is exactly singular.  in */
139/*               this case, info is the index of the first zero */
140/*               diagonal element of r and b is left unaltered. */
141
142/*     the parameters qy, qty, b, rsd, and xb are not referenced */
143/*     if their computation is not requested and in this case */
144/*     can be replaced by dummy variables in the calling program. */
145/*     to save storage, the user may in some cases use the same */
146/*     array for different parameters in the calling sequence.  a */
147/*     frequently occuring example is when one wishes to compute */
148/*     any of b, rsd, or xb and does not need y or qty.  in this */
149/*     case one may identify y, qty, and one of b, rsd, or xb, while */
150/*     providing separate arrays for anything else that is to be */
151/*     computed.  thus the calling sequence */
152
153/*          call dqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info) */
154
155/*     will result in the computation of b and rsd, with rsd */
156/*     overwriting y.  more generally, each item in the following */
157/*     list contains groups of permissible identifications for */
158/*     a single callinng sequence. */
159
160/*          1. (y,qty,b) (rsd) (xb) (qy) */
161
162/*          2. (y,qty,rsd) (b) (xb) (qy) */
163
164/*          3. (y,qty,xb) (b) (rsd) (qy) */
165
166/*          4. (y,qy) (qty,b) (rsd) (xb) */
167
168/*          5. (y,qy) (qty,rsd) (b) (xb) */
169
170/*          6. (y,qy) (qty,xb) (b) (rsd) */
171
172/*     in any group the value returned in the array allocated to */
173/*     the group corresponds to the last member of the group. */
174
175/*     linpack. this version dated 08/14/78 . */
176/*     g.w. stewart, university of maryland, argonne national lab. */
177
178/*     dqrsl uses the following functions and subprograms. */
179
180/*     blas daxpy,dcopy,ddot */
181/*     fortran dabs,min0,mod */
182
183/*     internal variables */
184
185
186
187/*     set info flag. */
188
189    /* Parameter adjustments */
190    x_dim1 = *ldx;
191    x_offset = 1 + x_dim1;
192    x -= x_offset;
193    --qraux;
194    --y;
195    --qy;
196    --qty;
197    --b;
198    --rsd;
199    --xb;
200
201    /* Function Body */
202    *info = 0;
203
204/*     determine what is to be computed. */
205
206    cqy = *job / 10000 != 0;
207    cqty = *job % 10000 != 0;
208    cb = *job % 1000 / 100 != 0;
209    cr = *job % 100 / 10 != 0;
210    cxb = *job % 10 != 0;
211/* Computing MIN */
212    i__1 = *k, i__2 = *n - 1;
213    ju = min(i__1,i__2);
214
215/*     special action when n=1. */
216
217    if (ju != 0) {
218    goto L40;
219    }
220    if (cqy) {
221    qy[1] = y[1];
222    }
223    if (cqty) {
224    qty[1] = y[1];
225    }
226    if (cxb) {
227    xb[1] = y[1];
228    }
229    if (! cb) {
230    goto L30;
231    }
232    if (x[x_dim1 + 1] != 0.) {
233    goto L10;
234    }
235    *info = 1;
236    goto L20;
237L10:
238    b[1] = y[1] / x[x_dim1 + 1];
239L20:
240L30:
241    if (cr) {
242    rsd[1] = 0.;
243    }
244    goto L250;
245L40:
246
247/*        set up to compute qy or qty. */
248
249    if (cqy) {
250    dcopy_(n, &y[1], &c__1, &qy[1], &c__1);
251    }
252    if (cqty) {
253    dcopy_(n, &y[1], &c__1, &qty[1], &c__1);
254    }
255    if (! cqy) {
256    goto L70;
257    }
258
259/*           compute qy. */
260
261    i__1 = ju;
262    for (jj = 1; jj <= i__1; ++jj) {
263    j = ju - jj + 1;
264    if (qraux[j] == 0.) {
265        goto L50;
266    }
267    temp = x[j + j * x_dim1];
268    x[j + j * x_dim1] = qraux[j];
269    i__2 = *n - j + 1;
270    t = -ddot_(&i__2, &x[j + j * x_dim1], &c__1, &qy[j], &c__1) / x[j + j
271        * x_dim1];
272    i__2 = *n - j + 1;
273    daxpy_(&i__2, &t, &x[j + j * x_dim1], &c__1, &qy[j], &c__1);
274    x[j + j * x_dim1] = temp;
275L50:
276/* L60: */
277    ;
278    }
279L70:
280    if (! cqty) {
281    goto L100;
282    }
283
284/*           compute trans(q)*y. */
285
286    i__1 = ju;
287    for (j = 1; j <= i__1; ++j) {
288    if (qraux[j] == 0.) {
289        goto L80;
290    }
291    temp = x[j + j * x_dim1];
292    x[j + j * x_dim1] = qraux[j];
293    i__2 = *n - j + 1;
294    t = -ddot_(&i__2, &x[j + j * x_dim1], &c__1, &qty[j], &c__1) / x[j +
295        j * x_dim1];
296    i__2 = *n - j + 1;
297    daxpy_(&i__2, &t, &x[j + j * x_dim1], &c__1, &qty[j], &c__1);
298    x[j + j * x_dim1] = temp;
299L80:
300/* L90: */
301    ;
302    }
303L100:
304
305/*        set up to compute b, rsd, or xb. */
306
307    if (cb) {
308    dcopy_(k, &qty[1], &c__1, &b[1], &c__1);
309    }
310    kp1 = *k + 1;
311    if (cxb) {
312    dcopy_(k, &qty[1], &c__1, &xb[1], &c__1);
313    }
314    if (cr && *k < *n) {
315    i__1 = *n - *k;
316    dcopy_(&i__1, &qty[kp1], &c__1, &rsd[kp1], &c__1);
317    }
318    if (! cxb || kp1 > *n) {
319    goto L120;
320    }
321    i__1 = *n;
322    for (i__ = kp1; i__ <= i__1; ++i__) {
323    xb[i__] = 0.;
324/* L110: */
325    }
326L120:
327    if (! cr) {
328    goto L140;
329    }
330    i__1 = *k;
331    for (i__ = 1; i__ <= i__1; ++i__) {
332    rsd[i__] = 0.;
333/* L130: */
334    }
335L140:
336    if (! cb) {
337    goto L190;
338    }
339
340/*           compute b. */
341
342    i__1 = *k;
343    for (jj = 1; jj <= i__1; ++jj) {
344    j = *k - jj + 1;
345    if (x[j + j * x_dim1] != 0.) {
346        goto L150;
347    }
348    *info = j;
349/*           ......exit */
350    goto L180;
351L150:
352    b[j] /= x[j + j * x_dim1];
353    if (j == 1) {
354        goto L160;
355    }
356    t = -b[j];
357    i__2 = j - 1;
358    daxpy_(&i__2, &t, &x[j * x_dim1 + 1], &c__1, &b[1], &c__1);
359L160:
360/* L170: */
361    ;
362    }
363L180:
364L190:
365    if (! cr && ! cxb) {
366    goto L240;
367    }
368
369/*           compute rsd or xb as required. */
370
371    i__1 = ju;
372    for (jj = 1; jj <= i__1; ++jj) {
373    j = ju - jj + 1;
374    if (qraux[j] == 0.) {
375        goto L220;
376    }
377    temp = x[j + j * x_dim1];
378    x[j + j * x_dim1] = qraux[j];
379    if (! cr) {
380        goto L200;
381    }
382    i__2 = *n - j + 1;
383    t = -ddot_(&i__2, &x[j + j * x_dim1], &c__1, &rsd[j], &c__1) / x[j +
384        j * x_dim1];
385    i__2 = *n - j + 1;
386    daxpy_(&i__2, &t, &x[j + j * x_dim1], &c__1, &rsd[j], &c__1);
387L200:
388    if (! cxb) {
389        goto L210;
390    }
391    i__2 = *n - j + 1;
392    t = -ddot_(&i__2, &x[j + j * x_dim1], &c__1, &xb[j], &c__1) / x[j + j
393        * x_dim1];
394    i__2 = *n - j + 1;
395    daxpy_(&i__2, &t, &x[j + j * x_dim1], &c__1, &xb[j], &c__1);
396L210:
397    x[j + j * x_dim1] = temp;
398L220:
399/* L230: */
400    ;
401    }
402L240:
403L250:
404    return 0;
405} /* dqrsl_ */
406
Note: See TracBrowser for help on using the repository browser.