/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:41 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dhfti.h" #include #include /* PARAMETER translations */ #define COL TRUE #define FACTOR 1000.0e0 #define ONE 1.0e0 #define ROW FALSE #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dhfti( double *a, long lda, long m1, long n1, double *b, long ldb, long kb, double tau, long *krank, double rnorm[], double work[], long ip[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) #define B(I_,J_) (*(b+(I_)*(ldb)+(J_))) LOGICAL32 comsqr; long int _l0, i, ii, j, jb, k, kp1, l, ldiag, lmax_, m, n, nterms; double hfac, sm1, small, tmp, uparam; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Ip = &ip[0] - 1; double *const Rnorm = &rnorm[0] - 1; double *const Work = &work[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2006-03-20 DHFTI Krogh Changed LDB to LDA in an error message. *>> 1996-03-30 DHFTI Krogh Added external statement. *>> 1994-10-20 DHFTI Krogh Changes to use M77CON *>> 1994-04-20 DHFTI CLL Edited to make DP & SP files similar. *>> 1993-02-09 CLL. Fixed index in 1st reference to [D/S]NRM2. *>> 1992-03-13 DHFTI FTK Removed implicit statements. *>> 1987-11-24 DHFTI Lawson Initial code. *--D replaces "?": ?HFTI, ?HTCC, ?HTGEN, ?DOT, ?NRM2 * * ------------------------------------------------------------------ * This subr solves the least squares problem * * A * X ~=~ B * * where A is a given M x N matrix, B is a given M x KB matrix and * X is the N x KB solution matrix to be determined. This includes * the usual special case of KB = 1 where B is an M-vector and the * solution, X, is an N-vector. * * This subr permits M > N, M = N, or M < N. This subr * determines the "pseudorank", i.e. the estimated rank, of A based * on a user-provided tolerance. If the pseudorank is less than N, * the minimal length solution, i.e. the pseudoinverse solution, to * the problem is computed. * * Note that this subr can be used to compute the pseudoinverse * of a matrix, A. Set B to the M x M identity matrix and the * solution matrix, X, will be the pseudoinverse of A. * * The algorithm is HFTI from the L & H book. This method does * a Householder QR decomposition from the left. Then if the * pseudorank is less than N it does a second Householder QR * decomposition from the right. * * The results returned in A(,), RNORM(), and IP() can be used * by subroutine SCOV1 or DCOV1 to compute the covariance matrix of * the solution vectors. * ------------------------------------------------------------------ * SUBROUTINE ARGUMENTS * * A(,) (In/Out) On input, contains the M x N matrix, A. Permit * M > N, M = N, or M < N. On return A(,) will contain an * upper triangular matrix of order KRANK that can be used * by subr _COV2 to compute a covariance matrix when * KRANK = N. * * LDA (In) The first dimension of the array A(,). * Require LDA .ge. M. * * M (In) No. of rows of matrices A and B. Require M .ge. 1. * * N (In) No. of columns of matrix A, and rows of matrix X. * Require N .ge. 0. * * B(,) (In/Out) If KB > 0, the array B(,) must initially * contain the right-side matrix, B, having M rows and KB * columns. On return the array B(,) will contain the * N x KB solution matrix X. * If KB = 0, this subr will not reference the array B(,). * * LDB (In) First dimensioning parameter for the array B(,). * If KB > 0, require LDB .ge. Max( M, N). * If KB = 0, require LDB .ge. 1. * * KB (In) No. of columns of the matrices B and X. * Require KB .ge. 0. * If KB = 0, this subr will not reference the array B(,). * * TAU (In) Absolute tolerance parameter provided by user for * pseudorank determination. * * KRANK (Out) Set by subr to indicate the pseudorank of A. * This means that the first KRANK diagonal elements in the * the upper triangular factor matrix derived from A each * exceed TAU in magnitude. Either KRANK = Min( M, N), or * the the magnitude of the diagonal element in position * KRANK + 1 is less than or equal to TAU. * * RNORM() (Out) On return, RNORM(J) will contain the euclidean * norm of the residual vector for the problem defined by * the Jth column vector of the input matrix, B, for * J = 1, ..., KB. * * WORK() (Scratch) Array used for work space by this subr. * Must be of length at least N. * * IP() (Work/Out) Integer array of length at least N in which * the subr will store column permutation information. * ----------------------------------------------------------------- * Subprograms referenced directly: ERMSG, ERMOR, IERM1, IERV1 * D1MACH, DHTCC, DHTGEN, DDOT, DNRM2 * Other subprograms needed: ERFIN * ----------------------------------------------------------------- * This code was originally developed by Charles L. Lawson and * Richard J. Hanson at Jet Propulsion Laboratory in 1973. The * original code was described and listed in the book, * */ /* Solving Least Squares Problems * C. L. Lawson and R. J. Hanson * Prentice-Hall, 1974 * * Feb 1985, Mar 1987, June 1987. C. L. Lawson & S. Y. Chiu, JPL. * Adapted code from the Lawson & Hanson book to Fortran 77 for use * in the JPL MATH77 library. * Changed code to provide oveflow avoidance. * Replaced previous scratch arrays H() and G() by WORK(). * Prefixing subprogram names with S or D for s.p. or d.p. versions. * Using generic names for intrinsic functions. * Adding calls to BLAS and MATH77 error processing subrs in some * program units. * ------------------------------------------------------------------ * 1983 Sept 22. CLL added computation of RNORM() for the * exceptional case of N = 0. * ----------------------------------------------------------------- */ /* ----------------------------------------------------------------- * */ m = m1; n = n1; if (((m < 1 || n < 0) || kb < 0) || lda < m) { ermsg( "DHFTI", 1, 0, "Bad argument values. Require M .ge. 1, N .ge. 0," , ',' ); ermor( "KB .ge. 0, and LDA .ge. M", ',' ); ierv1( "M", m, ',' ); ierv1( "N", n, ',' ); ierv1( "KB", kb, ',' ); ierv1( "LDA", lda, '.' ); *krank = 0; return; } else if (kb == 0) { if (ldb <= 0) { ierm1( "DHFTI", 2, 0, "Require LDB .ge. 1 when KB .eq. 0" , "LDB", ldb, '.' ); *krank = 0; return; } } else if (ldb < max( m, n )) { ierm1( "DHFTI", 3, 0, "Require LDB .ge. max(M,N) when KB .ge. 1" , "KB", kb, ',' ); ierv1( "LDB", ldb, '.' ); *krank = 0; return; } if (n == 0) { for (j = 1; j <= kb; j++) { Rnorm[j] = dnrm2( m, &B(j - 1,0), 1 ); } *krank = 0; return; } /* Here we have M > 0 and N > 0. */ small = FACTOR*DBL_EPSILON; ldiag = min( m, n ); for (j = 1; j <= ldiag; j++) { if (j == n) { /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Special for J = N. This case is much simpler than J < N * since there are no more columns of A beyond the jth to be * considered for interchange or to be triangularized. * */ Ip[n] = n; dhtcc( 1, n, n + 1, m, &A(n - 1,0), &uparam, b, ldb, kb ); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ } else { /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Here we have J < N. */ if (j == 1) { comsqr = TRUE; } else { /* Update scaled squared column lengths and set LMAX. * */ lmax_ = j; for (l = j; l <= n; l++) { Work[l] += -powi(hfac*A(l - 1,j - 2),2); if (Work[l] > Work[lmax_]) lmax_ = l; } comsqr = Work[lmax_] <= small; } if (comsqr) { /* Compute scaled squared column lengths and set LMAX. * Scaling using HFAC protects against overflow of squared * numbers. * */ nterms = m - j + 1; lmax_ = j; for (l = j; l <= n; l++) { Work[l] = dnrm2( nterms, &A(l - 1,j - 1), 1 ); if (Work[l] > Work[lmax_]) lmax_ = l; } if (Work[lmax_] == ZERO) { hfac = ONE; } else { hfac = ONE/Work[lmax_]; } for (l = j; l <= n; l++) { Work[l] = powi(hfac*Work[l],2); } } /* DO COLUMN INTERCHANGES IF NEEDED. * */ Ip[j] = lmax_; if (Ip[j] != j) { for (i = 1; i <= m; i++) { tmp = A(j - 1,i - 1); A(j - 1,i - 1) = A(lmax_ - 1,i - 1); A(lmax_ - 1,i - 1) = tmp; } Work[lmax_] = Work[j]; } /* Compute the J-th transformation and apply it to A and B. * Since we treated J = N as a special case we here have J < N * so the reference to A(1,J+1) is valid. * */ dhtcc( 1, j, j + 1, m, &A(j - 1,0), &uparam, &A(j,0), lda, n - j ); dhtcc( 2, j, j + 1, m, &A(j - 1,0), &uparam, b, ldb, kb ); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ } } /* DETERMINE THE PSEUDORANK, K, USING THE TOLERANCE, TAU. * */ k = ldiag; for (j = 1; j <= ldiag; j++) { if (fabs( A(j - 1,j - 1) ) <= tau) { k = j - 1; goto L_100; } } L_100: ; kp1 = k + 1; /* COMPUTE THE NORMS OF THE RESIDUAL VECTORS. * */ for (jb = 1; jb <= kb; jb++) { tmp = ZERO; for (i = kp1; i <= m; i++) { tmp += SQ(B(jb - 1,i - 1)); } Rnorm[jb] = sqrt( tmp ); } /* Special termination when Pseudorank = 0 */ if (k == 0) { for (jb = 1; jb <= kb; jb++) { for (i = 1; i <= n; i++) { B(jb - 1,i - 1) = ZERO; } } *krank = 0; return; } /* IF THE PSEUDORANK IS LESS THAN N COMPUTE HOUSEHOLDER * DECOMPOSITION OF FIRST K ROWS. * */ if (k != n) { for (ii = 1; ii <= k; ii++) { i = kp1 - ii; dhtgen( 1, i, kp1, n, &A(0,i - 1), lda, ROW, &Work[i], a, lda, i - 1, ROW ); } } for (jb = 1; jb <= kb; jb++) { /* SOLVE THE K BY K TRIANGULAR SYSTEM. * */ for (l = 1; l <= k; l++) { i = kp1 - l; if (i < k) { sm1 = ddot( k - i, &A(i,i - 1), lda, &B(jb - 1,i), 1 ); } else { sm1 = ZERO; } B(jb - 1,i - 1) = (B(jb - 1,i - 1) - sm1)/A(i - 1,i - 1); } /* COMPLETE COMPUTATION OF SOLUTION VECTOR. * .. */ if (k != n) { for (j = kp1; j <= n; j++) { B(jb - 1,j - 1) = ZERO; } for (i = 1; i <= k; i++) { dhtgen( 2, i, kp1, n, &A(0,i - 1), lda, ROW, &Work[i], &B(jb - 1,0), ldb, 1, COL ); } } /* RE-ORDER THE SOLUTION VECTOR TO COMPENSATE FOR THE * COLUMN INTERCHANGES. * */ for (j = ldiag; j >= 1; j--) { if (Ip[j] != j) { l = Ip[j]; tmp = B(jb - 1,l - 1); B(jb - 1,l - 1) = B(jb - 1,j - 1); B(jb - 1,j - 1) = tmp; } } } *krank = k; return; #undef B #undef A } /* end of function */