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

Revision 8069:5e19a5dcaac0, 8.3 KB checked in by ales_erjavec <ales.erjavec@…>, 3 years ago (diff)
Line 
1/* dqrdc2.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 "f2c.h"
14#include "linpack.h"
15
16/* Table of constant values */
17
18static integer c__1 = 1;
19
20
21/*     dqrdc2 uses householder transformations to compute the qr */
22/*     factorization of an n by p matrix x.  a limited column */
23/*     pivoting strategy based on the 2-norms of the reduced columns */
24/*     moves columns with near-zero norm to the right-hand edge of */
25/*     the x matrix.  this strategy means that sequential one */
26/*     degree-of-freedom effects can be computed in a natural way. */
27
28/*     i am very nervous about modifying linpack code in this way. */
29/*     if you are a computational linear algebra guru and you really */
30/*     understand how to solve this problem please feel free to */
31/*     suggest improvements to this code. */
32
33/*     Another change was to compute the rank. */
34
35/*     on entry */
36
37/*        x       double precision(ldx,p), where ldx .ge. n. */
38/*                x contains the matrix whose decomposition is to be */
39/*                computed. */
40
41/*        ldx     integer. */
42/*                ldx is the leading dimension of the array x. */
43
44/*        n       integer. */
45/*                n is the number of rows of the matrix x. */
46
47/*        p       integer. */
48/*                p is the number of columns of the matrix x. */
49
50/*        tol     double precision */
51/*                tol is the nonnegative tolerance used to */
52/*                determine the subset of the columns of x */
53/*                included in the solution. */
54
55/*        jpvt    integer(p). */
56/*                integers which are swapped in the same way as the */
57/*                the columns of x during pivoting.  on entry these */
58/*                should be set equal to the column indices of the */
59/*                columns of the x matrix (typically 1 to p). */
60
61/*        work    double precision(p,2). */
62/*                work is a work array. */
63
64/*     on return */
65
66/*        x       x contains in its upper triangle the upper */
67/*                triangular matrix r of the qr factorization. */
68/*                below its diagonal x contains information from */
69/*                which the orthogonal part of the decomposition */
70/*                can be recovered.  note that if pivoting has */
71/*                been requested, the decomposition is not that */
72/*                of the original matrix x but that of x */
73/*                with its columns permuted as described by jpvt. */
74
75/*        k       integer. */
76/*                k contains the number of columns of x judged */
77/*                to be linearly independent. */
78
79/*        qraux   double precision(p). */
80/*                qraux contains further information required to recover */
81/*                the orthogonal part of the decomposition. */
82
83/*        jpvt    jpvt(k) contains the index of the column of the */
84/*                original matrix that has been interchanged into */
85/*                the k-th column. */
86
87/*     original (dqrdc.f) linpack version dated 08/14/78 . */
88/*     g.w. stewart, university of maryland, argonne national lab. */
89
90/*     this version dated 22 august 1995 */
91/*     ross ihaka */
92
93/*     bug fixes 29 September 1999 BDR (p > n case, inaccurate ranks) */
94
95
96/*     dqrdc2 uses the following functions and subprograms. */
97
98/*     blas daxpy,ddot,dscal,dnrm2 */
99/*     fortran dabs,dmax1,min0,dsqrt */
100
101/* Subroutine */ int dqrdc2_(doublereal *x, integer *ldx, integer *n, integer
102    *p, doublereal *tol, integer *k, doublereal *qraux, integer *jpvt, 
103    doublereal *work)
104{
105    /* System generated locals */
106    integer x_dim1, x_offset, work_dim1, work_offset, i__1, i__2, i__3;
107    doublereal d__1, d__2;
108
109    /* Builtin functions */
110    double d_sign(doublereal *, doublereal *), sqrt(doublereal);
111
112    /* Local variables */
113    static integer i__, j, l;
114    static doublereal t, tt;
115    static integer lp1, lup;
116    static doublereal ttt;
117    extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
118        integer *), dnrm2_(integer *, doublereal *, integer *);
119    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
120        integer *), daxpy_(integer *, doublereal *, doublereal *, integer
121        *, doublereal *, integer *);
122    static doublereal nrmxl;
123
124
125/*     internal variables */
126
127
128
129/*     compute the norms of the columns of x. */
130
131    /* Parameter adjustments */
132    x_dim1 = *ldx;
133    x_offset = 1 + x_dim1;
134    x -= x_offset;
135    work_dim1 = *p;
136    work_offset = 1 + work_dim1;
137    work -= work_offset;
138    --qraux;
139    --jpvt;
140
141    /* Function Body */
142    i__1 = *p;
143    for (j = 1; j <= i__1; ++j) {
144    qraux[j] = dnrm2_(n, &x[j * x_dim1 + 1], &c__1);
145    work[j + work_dim1] = qraux[j];
146    work[j + (work_dim1 << 1)] = qraux[j];
147    if (work[j + (work_dim1 << 1)] == 0.) {
148        work[j + (work_dim1 << 1)] = 1.;
149    }
150/* L70: */
151    }
152
153/*     perform the householder reduction of x. */
154
155    lup = min(*n,*p);
156    *k = *p + 1;
157    i__1 = lup;
158    for (l = 1; l <= i__1; ++l) {
159
160/*     previous version only cycled l to lup */
161
162/*     cycle the columns from l to p left-to-right until one */
163/*     with non-negligible norm is located.  a column is considered */
164/*     to have become negligible if its norm has fallen below */
165/*     tol times its original norm.  the check for l .le. k */
166/*     avoids infinite cycling. */
167
168L80:
169    if (l >= *k || qraux[l] >= work[l + (work_dim1 << 1)] * *tol) {
170        goto L120;
171    }
172    lp1 = l + 1;
173    i__2 = *n;
174    for (i__ = 1; i__ <= i__2; ++i__) {
175        t = x[i__ + l * x_dim1];
176        i__3 = *p;
177        for (j = lp1; j <= i__3; ++j) {
178        x[i__ + (j - 1) * x_dim1] = x[i__ + j * x_dim1];
179/* L90: */
180        }
181        x[i__ + *p * x_dim1] = t;
182/* L100: */
183    }
184    i__ = jpvt[l];
185    t = qraux[l];
186    tt = work[l + work_dim1];
187    ttt = work[l + (work_dim1 << 1)];
188    i__2 = *p;
189    for (j = lp1; j <= i__2; ++j) {
190        jpvt[j - 1] = jpvt[j];
191        qraux[j - 1] = qraux[j];
192        work[j - 1 + work_dim1] = work[j + work_dim1];
193        work[j - 1 + (work_dim1 << 1)] = work[j + (work_dim1 << 1)];
194/* L110: */
195    }
196    jpvt[*p] = i__;
197    qraux[*p] = t;
198    work[*p + work_dim1] = tt;
199    work[*p + (work_dim1 << 1)] = ttt;
200    --(*k);
201    goto L80;
202L120:
203    if (l == *n) {
204        goto L190;
205    }
206
207/*           compute the householder transformation for column l. */
208
209    i__2 = *n - l + 1;
210    nrmxl = dnrm2_(&i__2, &x[l + l * x_dim1], &c__1);
211    if (nrmxl == 0.) {
212        goto L180;
213    }
214    if (x[l + l * x_dim1] != 0.) {
215        nrmxl = d_sign(&nrmxl, &x[l + l * x_dim1]);
216    }
217    i__2 = *n - l + 1;
218    d__1 = 1. / nrmxl;
219    dscal_(&i__2, &d__1, &x[l + l * x_dim1], &c__1);
220    x[l + l * x_dim1] += 1.;
221
222/*              apply the transformation to the remaining columns, */
223/*              updating the norms. */
224
225    lp1 = l + 1;
226    if (*p < lp1) {
227        goto L170;
228    }
229    i__2 = *p;
230    for (j = lp1; j <= i__2; ++j) {
231        i__3 = *n - l + 1;
232        t = -ddot_(&i__3, &x[l + l * x_dim1], &c__1, &x[l + j * x_dim1], &
233            c__1) / x[l + l * x_dim1];
234        i__3 = *n - l + 1;
235        daxpy_(&i__3, &t, &x[l + l * x_dim1], &c__1, &x[l + j * x_dim1], &
236            c__1);
237        if (qraux[j] == 0.) {
238        goto L150;
239        }
240/* Computing 2nd power */
241        d__2 = (d__1 = x[l + j * x_dim1], abs(d__1)) / qraux[j];
242        tt = 1. - d__2 * d__2;
243        tt = max(tt,0.);
244        t = tt;
245
246/* modified 9/99 by BDR. Re-compute norms if there is large reduction */
247/* The tolerance here is on the squared norm */
248/* In this version we need accurate norms, so re-compute often. */
249/*  work(j,1) is only updated in one case: looks like a bug -- no longer used */
250
251/*                     tt = 1.0d0 + 0.05d0*tt*(qraux(j)/work(j,1))**2 */
252/*                     if (tt .eq. 1.0d0) go to 130 */
253        if (abs(t) < 1e-6) {
254        goto L130;
255        }
256        qraux[j] *= sqrt(t);
257        goto L140;
258L130:
259        i__3 = *n - l;
260        qraux[j] = dnrm2_(&i__3, &x[l + 1 + j * x_dim1], &c__1);
261        work[j + work_dim1] = qraux[j];
262L140:
263L150:
264/* L160: */
265        ;
266    }
267L170:
268
269/*              save the transformation. */
270
271    qraux[l] = x[l + l * x_dim1];
272    x[l + l * x_dim1] = -nrmxl;
273L180:
274L190:
275/* L200: */
276    ;
277    }
278/* Computing min */
279    i__1 = *k - 1;
280    *k = min(i__1,*n);
281    return 0;
282} /* dqrdc2_ */
283
Note: See TracBrowser for help on using the repository browser.