/* Last edited on 2011-05-15 01:12:11 by stolfi */ /* Last edited on DATE TIME by USER */ /* qrfac.f -- translated by f2c (version of 17 January 1992 0:17:58). You must link the resulting object file with the libraries: -lf77 -li77 -lm -lc (in that order) */ #include #include #include /* Table of constant values */ static int c__1 = 1; /* Subroutine */ int qrfac_(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm, wa) int *m, *n; double *a; int *lda; logical *pivot; int *ipvt, *lipvt; double *rdiag, *acnorm, *wa; { /* Initialized data */ static double one = 1.; static double p05 = .05; static double zero = 0.; /* System generated locals */ int a_dim1, a_offset, i__1, i__2, i__3; double d__1, d__2, d__3; /* Builtin functions */ double sqrt(); /* Local variables */ static int kmax; static double temp; static int i, j, k, minmn; extern double enorm_(); static double epsmch; extern double dpmpar_(); static double ajnorm; static int jp1; static double sum; /* ********** */ /* subroutine qrfac */ /* this subroutine uses householder transformations with column */ /* pivoting (optional) to compute a qr factorization of the */ /* m by n matrix a. that is, qrfac determines an orthogonal */ /* matrix q, a permutation matrix p, and an upper trapezoidal */ /* matrix r with diagonal elements of nonincreasing magnitude, */ /* such that a*p = q*r. the householder transformation for */ /* column k, k = 1,2,...,min(m,n), is of the form */ /* t */ /* i - (1/u(k))*u*u */ /* where u has zeros in the first k-1 positions. the form of */ /* this transformation and the method of pivoting first */ /* appeared in the corresponding linpack subroutine. */ /* the subroutine statement is */ /* subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa) */ /* where */ /* m is a positive int input variable set to the number */ /* of rows of a. */ /* n is a positive int input variable set to the number */ /* of columns of a. */ /* a is an m by n array. on input a contains the matrix for */ /* which the qr factorization is to be computed. on output */ /* the strict upper trapezoidal part of a contains the strict */ /* upper trapezoidal part of r, and the lower trapezoidal */ /* part of a contains a factored form of q (the non-trivial */ /* elements of the u vectors described above). */ /* lda is a positive int input variable not less than m */ /* which specifies the leading dimension of the array a. */ /* pivot is a logical input variable. if pivot is set true, */ /* then column pivoting is enforced. if pivot is set false, */ /* then no column pivoting is done. */ /* ipvt is an int output array of length lipvt. ipvt */ /* defines the permutation matrix p such that a*p = q*r. */ /* column j of p is column ipvt(j) of the identity matrix. */ /* if pivot is false, ipvt is not referenced. */ /* lipvt is a positive int input variable. if pivot is false, */ /* then lipvt may be as small as 1. if pivot is true, then */ /* lipvt must be at least n. */ /* rdiag is an output array of length n which contains the */ /* diagonal elements of r. */ /* acnorm is an output array of length n which contains the */ /* norms of the corresponding columns of the input matrix a. */ /* if this information is not needed, then acnorm can coincide */ /* with rdiag. */ /* wa is a work array of length n. if pivot is false, then wa */ /* can coincide with rdiag. */ /* subprograms called */ /* minpack-supplied ... dpmpar,enorm */ /* fortran-supplied ... dmax1,dsqrt,min0 */ /* argonne national laboratory. minpack project. march 1980. */ /* burton s. garbow, kenneth e. hillstrom, jorge j. more */ /* ********** */ /* Parameter adjustments */ --wa; --acnorm; --rdiag; --ipvt; a_dim1 = *lda; a_offset = a_dim1 + 1; a -= a_offset; /* Function Body */ /* epsmch is the machine precision. */ epsmch = dpmpar_(&c__1); /* compute the initial column norms and initialize several arrays. */ i__1 = *n; for (j = 1; j <= i__1; ++j) { acnorm[j] = enorm_(m, &a[j * a_dim1 + 1]); rdiag[j] = acnorm[j]; wa[j] = rdiag[j]; if (*pivot) { ipvt[j] = j; } /* L10: */ } /* reduce a to r with householder transformations. */ minmn = min(*m,*n); i__1 = minmn; for (j = 1; j <= i__1; ++j) { if (! (*pivot)) { goto L40; } /* bring the column of largest norm into the pivot position. */ kmax = j; i__2 = *n; for (k = j; k <= i__2; ++k) { if (rdiag[k] > rdiag[kmax]) { kmax = k; } /* L20: */ } if (kmax == j) { goto L40; } i__2 = *m; for (i = 1; i <= i__2; ++i) { temp = a[i + j * a_dim1]; a[i + j * a_dim1] = a[i + kmax * a_dim1]; a[i + kmax * a_dim1] = temp; /* L30: */ } rdiag[kmax] = rdiag[j]; wa[kmax] = wa[j]; k = ipvt[j]; ipvt[j] = ipvt[kmax]; ipvt[kmax] = k; L40: /* compute the householder transformation to reduce the */ /* j-th column of a to a multiple of the j-th unit vector. */ i__2 = *m - j + 1; ajnorm = enorm_(&i__2, &a[j + j * a_dim1]); if (ajnorm == zero) { goto L100; } if (a[j + j * a_dim1] < zero) { ajnorm = -ajnorm; } i__2 = *m; for (i = j; i <= i__2; ++i) { a[i + j * a_dim1] /= ajnorm; /* L50: */ } a[j + j * a_dim1] += one; /* apply the transformation to the remaining columns */ /* and update the norms. */ jp1 = j + 1; if (*n < jp1) { goto L100; } i__2 = *n; for (k = jp1; k <= i__2; ++k) { sum = zero; i__3 = *m; for (i = j; i <= i__3; ++i) { sum += a[i + j * a_dim1] * a[i + k * a_dim1]; /* L60: */ } temp = sum / a[j + j * a_dim1]; i__3 = *m; for (i = j; i <= i__3; ++i) { a[i + k * a_dim1] -= temp * a[i + j * a_dim1]; /* L70: */ } if (! (*pivot) || rdiag[k] == zero) { goto L80; } temp = a[j + k * a_dim1] / rdiag[k]; /* Computing MAX */ /* Computing 2nd power */ d__3 = temp; d__1 = zero, d__2 = one - d__3 * d__3; rdiag[k] *= sqrt((max(d__1,d__2))); /* Computing 2nd power */ d__1 = rdiag[k] / wa[k]; if (p05 * (d__1 * d__1) > epsmch) { goto L80; } i__3 = *m - j; rdiag[k] = enorm_(&i__3, &a[jp1 + k * a_dim1]); wa[k] = rdiag[k]; L80: /* L90: */ ; } L100: rdiag[j] = -ajnorm; /* L110: */ } return 0; /* last card of subroutine qrfac. */ } /* qrfac_ */