/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:44 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dsvdrs.h" #include /* PARAMETER translations */ #define COL TRUE #define ONE 1.0e0 #define ROW FALSE #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dsvdrs( double *a, long lda, long m1, long n1, double *b, long ldb, long nb, double s[], double work[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) #define B(I_,J_) (*(b+(I_)*(ldb)+(J_))) long int i, ipass, j, k, l, lbase, mwork, np1, ns, nwork; double uparam; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const S = &s[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. *>> 2001-05-25 DSVDRS Krogh Minor change for making .f90 version. *>> 1999-08-12 DSVDRS Krogh Loop at 180 => 2 loops for Sun optimizer. *>> 1994-10-20 DSVDRS Krogh Changes to use M77CON *>> 1992-06-16 DSVDRS CLL *>> 1992-03-13 DSVDRS FTK Removed implicit statements. *>> 1989-10-20 CLL *>> 1989-06-05 DSVDRS Snyder Add END statement, polish a little. *>> 1987-11-24 DSVDRS Lawson Adaptation to Fortran 77 and MATH77. *>> 1974-01-01 Lawson and Hanson initial code. * * This subr computes the singular value decomposition of the * given M1 x N1 matrix, A, and optionally applys the transformations * from the left to the NB column vectors of the M1 x NB matrix B. * Either M1 .ge. N1 or M1 .lt. N1 is permitted. * * The singular value decomposition of A is of the form * * A = U * S * V**t * * where U is M1 x M1 orthogonal, S is M1 x N1 diagonal with the * diagonal terms nonnegative and ordered from large to small, and * V is N1 x N1 orthogonal. Note that these matrices also satisfy * * S = (U**t) * A * V * * The matrix V is returned in the leading N1 rows and * columns of the array A(,). * * The singular values, i.e. the diagonal terms of the matrix S, * are returned in the array S(). If M1 .lt. N1, positions M1+1 * through N1 of S() will be set to zero. * * The product matrix G = U**t * B replaces the given matrix B * in the array B(,). * * If the user wishes to obtain a minimum length least squares * solution of the linear system * * A * X ~=~ B * * the solution X can be constructed, following use of this subr, * by computing the sum for i = 1, ..., R of the outer products * * (Col i of V) * (1/S(i)) * (Row i of G) * * Here R denotes the pseudorank of A which the user may choose * in the range 0 through Min(M1, N1) based on the sizes of the * singular values. * ------------------------------------------------------------------ * Subroutine Arguments * * A(,) (In/Out) On input contains the M1 x N1 matrix A. * On output contains the N1 x N1 matrix V. * * LDA (In) First dimensioning parameter for A(,). * Require LDA .ge. Max(M1, N1). * * M1 (In) No. of rows of matrices A, B, and G. * Require M1 > 0. * * N1 (In) No. of cols of matrix A, No. of rows and cols of * matrix V. Permit M1 .ge. N1 or M1 .lt. N1. * Require N1 > 0. * * B(,) (In/Out) If NB .gt. 0 this array must contain an * M1 x NB matrix on input and will contain the * M1 x NB product matrix, G = (U**t) * B on output. * * LDB (In) First dimensioning parameter for B(,). * Require LDB .ge. M1. * * NB (In) No. of cols in the matrices B and G. * Require NB .ge. 0. * * S() (Out) Must be dimensioned at least N1. On return will * contain the singular values of A, with the ordering * S(1) .ge. S(2) .ge. ... .ge. S(N1) .ge. 0. * If M1 .lt. N1 the singular values indexed from M1+1 * through N1 will be zero. * If the given integer arguments are not consistent, this * subr will return immediately, setting S(1) = -1.0. * * WORK() (Scratch) Work space of total size at least 2*N1. * Locations 1 thru N1 will hold the off-diagonal terms of * the bidiagonal matrix for subr _QRBD. Locations N1+1 * thru 2*N1 will save info from one call to the next of * _HTGEN. * ------------------------------------------------------------------ * Subprograms referenced directly: DCOPY, DHTCC, DHTGEN,DQRBD,DSWAP, * ERMSG, ERMOR, IERV1 * Other subprograms needed: ERFIN, *1MACH (Fortran 77 only), *ROT, * *ROTG (* = S or D) * ------------------------------------------------------------------ * This subr gives special treatment to the cases of entire rows * and/or columns of A being zero. This provides the feature that * singular values that are zero because of the zero structure of A * will be produced as exact zeros rather than as relatively small */ /* nonzero values. This is convenient for cases in which a user * wishes to remove a variable from a problem by simply zeroing the * corresponding column of A. * * Method for this special feature: * * 1. EXCHANGE COLS OF A TO PACK NONZERO COLS TO THE LEFT. * SET NWORK = NO. OF NONZERO COLS. * USE LOCATIONS A(1,NWORK+1),...,A(1,N1) TO RECORD THE * COL PERMUTATIONS. * 2. EXCHANGE ROWS OF A TO PACK NONZERO ROWS TO THE TOP. * QUIT PACKING IF FIND NWORK NONZERO ROWS. MAKE SAME ROW * EXCHANGES IN B. SET MWORK SO THAT ALL NONZERO ROWS OF THE * PERMUTED A ARE IN FIRST MWORK ROWS. IF MWORK .LE. NWORK THEN * ALL MWORK ROWS ARE NONZERO. IF MWORK .GT. NWORK THEN THE * FIRST NWORK ROWS ARE KNOWN TO BE NONZERO, AND ROWS NWORK+1 * THRU MWORK MAY BE ZERO OR NONZERO. * 3. APPLY THE BASIC SVD ALGORITHM TO THE MWORK BY NWORK PROBLEM. * 4. MOVE PERMUTATION RECORD FROM A(,) TO S(I),I = NWORK+1,...,N1. * 5. BUILD V UP FROM NWORK BY NWORK TO N1 BY N1 BY PLACING ONES ON * THE DIAGONAL AND ZEROS ELSEWHERE. THIS IS ONLY PARTLY DONE * EXPLICITLY. IT IS COMPLETED DURING STEP 6. * 6. EXCHANGE ROWS OF V TO COMPENSATE FOR COL EXCHANGES OF STEP 2. * 7. PLACE ZEROS IN S(I),I = NWORK+1,...,N1 TO REPRESENT ZERO * SING VALS. * * ------------------------------------------------------------------ * 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, C. L. Lawson & S. Y. Chiu, JPL. Adapted code * from the Lawson & Hanson book to Fortran 77 for use in the * JPL MATH77 library. * 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. * 1989-10-20 CLL Moved integer declaration to earlier position to * avoid warning msg from Cray compiler. * 1992-06-16 CLL. Using "min" fcn in arg list of calls to _HTCC * and _HTGEN to avoid "index out of range" when using a bounds * checker, even though no reference would be made to the "out of * range" location. * ------------------------------------------------------------------ *--D replaces "?": ?SVDRS, ?COPY, ?HTCC, ?HTGEN, ?QRBD, ?SWAP * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ nwork = n1; lbase = nwork; if (((lda < max( m1, nwork ) || m1 <= 0) || nwork <= 0) || (nb > 0 && ldb < m1)) { ermsg( "DSVDRS", 1, 0, "Bad argument values. Require LDA .ge. max(M,N), M > 0," , ',' ); ermor( "N > 0, LDB .ge. M if NB > 0", ',' ); ierv1( "LDA", lda, ',' ); ierv1( "M", m1, ',' ); ierv1( "N", nwork, ',' ); ierv1( "LDB", ldb, ',' ); ierv1( "NB", nb, '.' ); S[1] = -ONE; return; } /* BEGIN.. SPECIAL FOR ZERO ROWS AND COLS. * * PACK THE NONZERO COLS TO THE LEFT * */ for (j = nwork; j >= 1; j--) { for (i = 1; i <= m1; i++) { if (A(j - 1,i - 1) != ZERO) goto L_50; } /* COL J IS ZERO. EXCHANGE IT WITH COL NWORK. * */ if (j != nwork) { for (i = 1; i <= m1; i++) { A(j - 1,i - 1) = A(nwork - 1,i - 1); } } A(nwork - 1,0) = j; nwork -= 1; L_50: ; } /* IF NWORK = 0 THEN A IS ENTIRELY ZERO AND SVD * COMPUTATION CAN BE SKIPPED */ ns = 0; if (nwork != 0) { /* --------------------------------------------------------------- * PACK NONZERO ROWS TO THE TOP * QUIT PACKING IF FIND NWORK NONZERO ROWS */ i = 1; mwork = m1; L_60: if (i <= nwork && i < mwork) { for (j = i; j <= nwork; j++) { if (A(j - 1,i - 1) != ZERO) goto L_145; } for (j = 1; j <= (i - 1); j++) { if (A(j - 1,i - 1) != ZERO) goto L_145; } /* ROW I IS ZERO * EXCHANGE ROWS I AND MWORK */ if (nb > 0) dswap( nb, &B(0,mwork - 1), ldb, &B(0,i - 1), ldb ); dcopy( nwork, &A(0,mwork - 1), lda, &A(0,i - 1), lda ); if (mwork <= nwork) { for (j = 1; j <= nwork; j++) { A(j - 1,mwork - 1) = ZERO; } } /* EXCHANGE IS FINISHED */ mwork -= 1; goto L_60; /* Here when row I is not all zero. No exchange needed. */ L_145: i += 1; goto L_60; } /* END.. SPECIAL FOR ZERO ROWS AND COLUMNS * --------------------------------------------------------------- * BEGIN.. SVD ALGORITHM * METHOD.. * (1) REDUCE THE MATRIX TO UPPER BIDIAGONAL FORM WITH * HOUSEHOLDER TRANSFORMATIONS. * H(NWORK)...H(1)AQ(1)...Q(NWORK-2) = (D**T,0)**T * WHERE D IS UPPER BIDIAGONAL. * * (2) APPLY H(NWORK)...H(1) TO B. HERE H(NWORK)...H(1)*B REPLAC * IN STORAGE. * * (3) THE MATRIX PRODUCT W = Q(1)...Q(NWORK-2) OVERWRITES THE F * NWORK ROWS OF A IN STORAGE. * * (4) AN SVD FOR D IS COMPUTED. HERE K ROTATIONS RI AND PI ARE * COMPUTED SO THAT * RK...R1*D*P1**(T)...PK**(T) = DIAG(S1,...,S.NS) * TO WORKING ACCURACY. THE SI ARE NONNEGATIVE AND NONINCREASING. * HERE RK...R1*B OVERWRITES B IN STORAGE WHILE * A*P1**(T)...PK**(T) OVERWRITES A IN STORAGE. * * (5) IT FOLLOWS THAT,WITH THE PROPER DEFINITIONS, * U**(T)*B OVERWRITES B, WHILE V OVERWRITES THE FIRST NWORK ROW AND * COLUMNS OF A. * ------------------------------------------------------------------ */ l = min( mwork, nwork ); /* THE FOLLOWING LOOP REDUCES A TO UPPER BIDIAGONAL AND * ALSO APPLIES THE PREMULTIPLYING TRANSFORMATIONS TO B. * */ for (j = 1; j <= l; j++) { if (j < mwork) { dhtcc( 1, j, j + 1, mwork, &A(j - 1,0), &uparam, &A(min( j + 1, nwork ) - 1,0), lda, nwork - j ); dhtcc( 2, j, j + 1, mwork, &A(j - 1,0), &uparam, b, ldb, nb ); } if (j < nwork - 1) { dhtgen( 1, j + 1, j + 2, nwork, &A(0,j - 1), lda, ROW, &Work[lbase + j], &A(0,min( j + 1, mwork ) - 1), lda, mwork - j, ROW ); } } /* Copy the bidiagonal matrix into the arrays * S() and WORK(1:N1) for _QRBD. */ for (j = 1; j <= l; j++) { S[j] = A(j - 1,j - 1); } for (j = 2; j <= l; j++) { Work[j] = A(j - 1,j - 2); } ns = nwork; if (mwork < nwork) { ns = mwork + 1; S[ns] = ZERO; Work[ns] = A(mwork,mwork - 1); } /* CONSTRUCT THE EXPLICIT NWORK BY NWORK PRODUCT MATRIX, * W = Q1*Q2*...*QL*I IN THE ARRAY A(). * */ for (i = nwork; i >= 1; i--) { if (i <= min( mwork, nwork - 2 )) dhtgen( 2, i + 1, i + 2, nwork, &A(0,i - 1), lda, ROW, &Work[lbase + i], &A(min( i + 1, nwork ) - 1,0), lda, nwork - i, COL ); for (j = 1; j <= nwork; j++) { A(j - 1,i - 1) = ZERO; } A(i - 1,i - 1) = ONE; } /* COMPUTE THE SVD OF THE BIDIAGONAL MATRIX * */ dqrbd( &ipass, s, &Work[1], ns, a, lda, nwork, b, ldb, nb ); /* If the above subr has a convergence failure it will report * the error via the MATH77 error message subrs and return * with IPASS = 2. Since the error has been reported, we don't * report it again here but just return. * */ if (ipass == 2) return; /* --------------------------------------------------------------- */ } for (j = ns + 1; j <= nwork; j++) { S[j] = ZERO; } if (nwork == n1) return; np1 = nwork + 1; /* MOVE RECORD OF PERMUTATIONS * AND STORE ZEROS */ for (j = np1; j <= n1; j++) { S[j] = A(j - 1,0); for (i = 1; i <= nwork; i++) { A(j - 1,i - 1) = ZERO; } } /* PERMUTE ROWS AND SET ZERO SINGULAR VALUES. */ for (k = np1; k <= n1; k++) { i = S[k]; S[k] = ZERO; for (j = 1; j <= n1; j++) { A(j - 1,k - 1) = A(j - 1,i - 1); A(j - 1,i - 1) = ZERO; } A(k - 1,i - 1) = ONE; } /* End.. Special for zero rows and columns. */ return; #undef B #undef A } /* end of function */