/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:57 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dsva.h" #include #include /* PARAMETER translations */ #define ONE 1.0e0 #define TEN 10.0e0 #define THOU 1000.0e0 #define TWENTY 20.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dsva( double *a, long mda, long m, long n, long mdata, double b[], double sing[], long kpvec[], char *names, int names_s, long iscale, double d[], double work[]) { #define A(I_,J_) (*(a+(I_)*(mda)+(J_))) #define NAMES(I_,J_) (names+(I_)*(names_s)+(J_)) LOGICAL32 blk[6], narrow; long int i, ie, ipass, j, k, minmn, minmn1, mpass, nsol, prblk, unit, width; double a1, a2, a3, a4, alamb, aln10, del, denom, el, el2, pcoef, rl, rnorm, rs, sb, sl, temp, yl, ynorm, ys, ysq; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const B = &b[0] - 1; LOGICAL32 *const Blk = &blk[0] - 1; double *const D = &d[0] - 1; long *const Kpvec = &kpvec[0] - 1; double *const Sing = &sing[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. *>> 1996-05-10 DSVA Krogh Only "*" output on C version. *>> 1995-11-15 DSVA Krogh Moved formats up for C conversion. *>> 1994-10-20 DSVA Krogh Changes to use M77CON *>> 1992-04-23 DSVA CLL Replaced Hollerith with '...' in a FORMAT. *>> 1992-03-18 CLL New treatment of WIDTH and length of labeling names. *>> 1989-11-21 CLL Fixed bug that caused divide by 0 when M < N. *>> 1989-03-07 DSVA CLL Changed logical PRINT to integer vector KPVEC. *>> 1988-12-01 DSVA CLL Added PRINT into arg list. *>> 1987-11-24 DSVA Lawson Initial code. * * Singular Value Analysis. Computes the singular value * decomposition of the matrix of a least squares problem, and * produces a printed report. * ------------------------------------------------------------------ * Subroutine Arguments * A(,) [inout] On entry, contains the M x N matrix of the least * squares problem to be analyzed. This could be a matrix * obtained by preliminary orthogonal transformations * applied to the actual problem matrix which may have had * more rows (See MDATA below.) * * MDA [in] First dimensioning parameter for A(,). Require * MDA .ge. max(M, N). * * M,N [in] No. of rows and columns, respectively, in the * matrix, A. Either M > N or M .le. N is permitted. * Require M > 0 and N > 0. * * MDATA [in] No. of rows in actual least squares problem. * Generally MDATA .ge. M. MDATA is used only in computing * statistics for the report and is not used as a loop * count or array dimension. * * B() [inout] On entry, contains the right-side vector, b, of the * least squares problem. This vector is of length, M. * On return, contains the vector, g = (U**t)*b, where U * comes from the singular value decomposition of A. The * vector , g, is also of length M. * * SING() [out] On return, contains the singular values of A, in * descending order, in locations indexed 1 thru min(M,N). * * KPVEC() [integer array, in] Array of integers to select print * options. KPVEC(1) determines whether the rest of * the array is to be used or ignored. * If KPVEC(1) = 1, the contents of (KPVEC(I), I=2,4) * will be used to set internal variables as follows: * PRBLK = KPVEC(2) * UNIT = KPVEC(3) * WIDTH = KPVEC(4) * If KPVEC(1) = 0 default settings will be used. The user * need not dimension KPVEC() greater than 1. The subr will * set PRBLK = 111111, UNIT = -1, and WIDTH = 69. * * The internal variables PRBLK, UNIT, and WIDTH are * interpreted as follows: * * PRBLK The decimal representation of PRBLK must be * representable as at most 6 digits, each being 0 or 1. * The decimal digits will be interpreted as independant * on/off flags for the 6 possible blocks of printed output. * Examples: 111111 selects all blocks, 0 suppresses all * printing, 101010 selects the 1st, 3rd, and 5th blocks, * etc. * The six blocks are: * 1. Header, with size and scaling option parameters. * 2. V-matrix. Amount of output depends on M and N. * 3. Singular values and related quantities. Amount of * output depends on N. * 4. Listing of YNORM and RNORM and their logarithms. * Amount of output depends on N. * 5. Levenberg-Marquart analysis. * 6. Candidate solutions. Amount of output depends on * M and N. * * UNIT Selects the output unit. If UNIT .ge. 0, * UNIT will be used as the output unit number. * If UNIT = -1, output will be written to the "*" output * unit, i.e., the standard system output unit. * The calling program unit is responsible for opening * and/or closing the selected output unit if the host * system requires these actions. * * WIDTH Default value is 79. Determines the width of * blocks 2, 3, and 6 of the output report. * Block 3 will use 95(+1) cols if WIDTH .ge. 95, and otherwise * 69(+1) cols. * Blocks 2 and 6 are printed by subroutine DPRTSV. These blocks * generally use at most WIDTH(+1) cols, but will use more if * the names are so long that more space is needed to print one * name and one numeric column. The (+1)'s above are reminders * that in all cases there is one extra initial column for Fortran * "carriage control". The carriage control character will always * be a blank. * Blocks 1, 4, and 5 have fixed widths of 63(+1), 66(+1) and * 66(+1), respectively. * */ /* NAMES() [in] NAMES(j), for j = 1, ..., N, may contain a * name for the jth component of the solution * vector. The declared length of the elements of the * NAMES() array is not specifically limited, but a * greater length reduces the space available for columns * of the matris to be printed. * If NAMES(1) contains only blank characters, * it will be assumed that no names have been provided, * and this subr will not access the NAMES() array beyond * the first element. * * ISCALE [in] Set by the user to 1, 2, or 3 to select the column * scaling option. * 1 SUBR WILL USE IDENTITY SCALING AND IGNORE THE D() * ARRAY. * 2 SUBR WILL SCALE NONZERO COLS TO HAVE UNIT EUCLID- * EAN LENGTH AND WILL STORE RECIPROCAL LENGTHS OF * ORIGINAL NONZERO COLS IN D(). * 3 USER SUPPLIES COL SCALE FACTORS IN D(). SUBR * WILL MULT COL J BY D(J) AND REMOVE THE SCALING * FROM THE SOLN AT THE END. * * D() [ignored or out or in] Usage of D() depends on ISCALE as * described above. When used, its length must be * at least N. * * WORK() [scratch] Work space of length at least 2*N. Used * directly in this subr and also in _SVDRS. * ------------------------------------------------------------------ *--D replaces "?": ?SVA, ?SVDRS, ?PRTSV * ------------------------------------------------------------------ * 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. * Dec 1988 CLL. For MATH77 this subr was intended to have PRINT in * the arg list, and was described that way in the MATH77 manual, but * PRINT was omitted in the code. Now fixing this at Rel. 2.2 by * adding PRINT to the arg list. * Dec 1988 CLL. Also for MATH77 the WORK() array was introduced. * Previous version used a long SING() array for work space. WORK() * was introduced incorrectly in the S.P. version. Now correcting * this and using the processor CHGTYP to keep the S.P. and D.P. * versions similar. * 1989-03-07 CLL> Changing the logical PRINT to an integer vector * KPVEC(). Giving user control of quantity, destination, and * line width of output. * ------------------------------------------------------------------ */ /*++ CODE for ~.C. is inactive * logical STAR *++ END * */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* ------------------------------------------------------------------ */ if (m <= 0 || n <= 0) return; minmn = min( m, n ); minmn1 = minmn + 1; if (Kpvec[1] == 0) { prblk = 111111; unit = -1; width = 79; } else { prblk = Kpvec[2]; unit = Kpvec[3]; width = Kpvec[4]; } /*++ CODE for ~.C. is inactive * STAR = UNIT .lt. 0 *++ END * Build logical array BLK() by testing * decimal digits of PRBLK. */ for (i = 6; i >= 1; i--) { j = prblk/10; Blk[i] = (prblk - 10*j) > 0; prblk = j; } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Optionally print header and M, N, MDATA */ if (Blk[1]) { /*++ CODE for ~.C. is inactive * if(STAR) then *++ END */ printf(" \n Singular Value Analysis of the least squares problem, A*X = B,\n scaled as (A*D)*Y = B.\n"); printf(" M = %6ld, N =%4ld, MDATA =%8ld\n", m, n, mdata); /*++ CODE for ~.C. is inactive * else * write (UNIT,270) * write (UNIT,260) M,N,MDATA * endif *++ END */ } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Handle scaling as selected by ISCALE. */ if (iscale == 1) { if (Blk[1]) { /*++ CODE for ~.C. is inactive * if(STAR) then *++ END */ printf(" \n Scaling option No. 1. D is the identity matrix.\n \n"); /*++ CODE for ~.C. is inactive * else * write (UNIT,290) * endif *++ END */ } } else { /* Apply column scaling to A. * */ for (j = 1; j <= n; j++) { a1 = D[j]; if (iscale <= 2) { sb = ZERO; for (i = 1; i <= m; i++) { sb += SQ(A(j - 1,i - 1)); } a1 = sqrt( sb ); if (a1 == ZERO) a1 = ONE; a1 = ONE/a1; D[j] = a1; } for (i = 1; i <= m; i++) { A(j - 1,i - 1) *= a1; } } if (Blk[1]) { /*++ CODE for ~.C. is inactive * if(STAR) then *++ END */ printf(" \n Scaling option No.%2ld", iscale); printf(". D is a diagonal matrix with the following diagonal elements..\n "); for (j = 1; j <= n; j++) { printf("%12.4e", D[j]); } printf("\n"); /*++ CODE for ~.C. is inactive * else * write (UNIT,280) ISCALE,(D(J),J = 1,N) * endif *++ END */ } } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Compute the Singular Value Decomposition of the scaled matrix. * */ dsvdrs( a, mda, m, n, b, m, 1, sing, work ); /* Determine NSOL. */ nsol = minmn; for (j = 1; j <= minmn; j++) { if (Sing[j] == ZERO) { nsol = j - 1; goto L_65; } } L_65: ; /* The array B() contains the vector G. * Compute cumulative sums of squares of components of * G and store them in WORK(I), I = 1,...,MINMN+1 * */ sb = ZERO; for (i = minmn1; i <= m; i++) { sb += SQ(B[i]); } Work[minmn + 1] = sb; for (j = minmn; j >= 1; j--) { sb += SQ(B[j]); Work[j] = sb; } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * PRINT THE V MATRIX. * */ if (Blk[2]) dprtsv( a, mda, n, n, names,names_s, 1, unit, width ); /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * REPLACE V BY D*V IN THE ARRAY A() */ if (iscale > 1) { for (i = 1; i <= n; i++) { for (j = 1; j <= n; j++) { A(j - 1,i - 1) *= D[i]; } } } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (Blk[3]) { /* Print singular values and other summary results. * * Output will be done using one of two layouts. The narrow * layout uses 69 cols + 1 for carriage control, and makes two passes * through the computation. * The wide layout uses 95 cols + 1 for carriage control, and makes * only one pass through the computation. * * G NOW IN B() ARRAY. V NOW IN A(,) ARRAY. * */ narrow = width < 95; mpass = 1; if (narrow) mpass = 2; for (ipass = 1; ipass <= mpass; ipass++) { /*++ CODE for ~.C. is inactive * if(STAR) then *++ END */ if (narrow) { if (ipass == 1) { printf(" \n INDEX SING. VAL. P COEF RECIPROCAL G COEF SCALED SQRT\n SING. VAL. of CUM.S.S.\n"); } else { printf(" \n INDEX SING. VAL. G COEF G**2 CUMULATIVE SCALED SQRT\n SUM of SQRS of CUM.S.S.\n"); } } else { printf(" \n INDEX SING. VAL. P COEF RECIPROCAL G COEF G**2 CUMULATIVE SCALED SQRT\n SING. VAL. SUM of SQRS of CUM.S.S.\n"); } /*++ CODE for ~.C. is inactive * else * if(NARROW) then * if(IPASS .eq. 1) then * write(UNIT,221) * else * write(UNIT,222) * endif * else * write (UNIT,220) * endif * endif *++ END * The following stmt converts from * integer to floating-point. */ denom = max( 1, mdata ); a3 = Work[1]; a4 = sqrt( a3/denom ); /*++ CODE for ~.C. is inactive * if(STAR) then *++ END */ if (narrow) { if (ipass == 1) { printf(" 0 %13.4g\n", a4); } else { printf(" 0 %13.4g%13.4g\n", a3, a4); } } else { printf(" 0 %13.4g%13.4g\n", a3, a4); } /*++ CODE for ~.C. is inactive * else * if(NARROW) then * if(IPASS .eq. 1) then * write(UNIT,231) A4 * else * write(UNIT,232) A3, A4 * endif * else * write (UNIT,230) A3,A4 * endif * endif *++ END * */ for (k = 1; k <= minmn; k++) { if (Sing[k] == ZERO) { pcoef = ZERO; /*++ CODE for ~.C. is inactive * if(STAR) then *++ END */ printf(" %5ld%12.4g\n", k, Sing[k]); /*++ CODE for ~.C. is inactive * else * write (UNIT,240) K,SING(K) * endif *++ END */ } else { pcoef = B[k]/Sing[k]; a1 = ONE/Sing[k]; a2 = SQ(B[k]); a3 = Work[k + 1]; /* The following stmt converts from * integer to floating-point. */ temp = max( 1, mdata - k ); a4 = sqrt( a3/temp ); /*++ CODE for ~.C. is inactive * if(STAR) then *++ END */ if (narrow) { if (ipass == 1) { printf(" %5ld%12.4g%13.4g%13.4g%13.4g%13.4g\n", k, Sing[k], pcoef, a1, B[k], a4); } else { printf(" %5ld%12.4g%13.4g%13.4g%13.4g%13.4g\n", k, Sing[k], B[k], a2, a3, a4); } } else { printf(" %5ld%12.4g%13.4g%13.4g%13.4g%13.4g%13.4g%13.4g\n", k, Sing[k], pcoef, a1, B[k], a2, a3, a4); } /*++ CODE for ~.C. is inactive * else * if(NARROW) then * if(IPASS .eq. 1) then * write(UNIT,240) K,SING(K),PCOEF,A1,B(K), A4 * else * write(UNIT,240) K,SING(K), B(K),A2,A3,A4 * endif * else * write (UNIT,240) K,SING(K),PCOEF,A1,B(K),A2,A3,A4 * endif */ /* endif *++ END */ } } } } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (Blk[4]) { /* Compute and print values of YNORM, RNORM and their logarithms. * *++ CODE for ~.C. is inactive * if(STAR) then *++ END */ printf(" \n INDEX YNORM RNORM LOG10 LOG10\n YNORM RNORM\n \n"); /*++ CODE for ~.C. is inactive * else * write (UNIT,300) * endif *++ END */ ysq = ZERO; for (j = 0; j <= nsol; j++) { if (j != 0) ysq += powi(B[j]/Sing[j],2); ynorm = sqrt( ysq ); rnorm = sqrt( Work[j + 1] ); yl = -THOU; if (ynorm > ZERO) yl = log10( ynorm ); rl = -THOU; if (rnorm > ZERO) rl = log10( rnorm ); /*++ CODE for ~.C. is inactive * if(STAR) then *++ END */ printf(" %5ld %11.3e%11.3e %11.3f%11.3f\n", j, ynorm, rnorm, yl, rl); /*++ CODE for ~.C. is inactive * else * write (UNIT,310) J,YNORM,RNORM,YL,RL * endif *++ END */ } } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (Blk[5] && Sing[1] != ZERO) { /* COMPUTE VALUES OF XNORM AND RNORM FOR A SEQUENCE OF VALUES OF * THE LEVENBERG-MARQUARDT PARAMETER. * */ el = log10( Sing[1] ) + ONE; el2 = log10( Sing[nsol] ) - ONE; del = (el2 - el)/TWENTY; aln10 = log( TEN ); /*++ CODE for ~.C. is inactive * if(STAR) then *++ END */ printf(" \n Norms of solution and residual vectors for a range of values\n of the Levenberg-Marquardt parameter, LAMBDA.\n\n LAMBDA YNORM RNORM LOG10 LOG10 LOG10\n LAMBDA YNORM RNORM\n"); /*++ CODE for ~.C. is inactive * else * write (UNIT,320) * endif *++ END */ for (ie = 1; ie <= 21; ie++) { /* COMPUTE ALAMB = 10.0**EL */ alamb = exp( aln10*el ); ys = ZERO; rs = Work[nsol + 1]; for (i = 1; i <= minmn; i++) { sl = SQ(Sing[i]) + SQ(alamb); ys += powi(B[i]*Sing[i]/sl,2); rs += powi(B[i]*(SQ(alamb))/sl,2); } ynorm = sqrt( ys ); rnorm = sqrt( rs ); rl = -THOU; if (rnorm > ZERO) rl = log10( rnorm ); yl = -THOU; if (ynorm > ZERO) yl = log10( ynorm ); /*++ CODE for ~.C. is inactive * if(STAR) then *++ END */ printf(" %11.3e%11.3e%11.3e%11.3f%11.3f%11.3f\n", alamb, ynorm, rnorm, el, yl, rl); /*++ CODE for ~.C. is inactive * else * write (UNIT,330) ALAMB,YNORM,RNORM,EL,YL,RL * endif *++ END */ el += del; } } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Compute and optionally print candidate solutions. * */ for (k = 1; k <= nsol; k++) { pcoef = B[k]/Sing[k]; for (i = 1; i <= n; i++) { A(k - 1,i - 1) *= pcoef; if (k > 1) A(k - 1,i - 1) += A(k - 2,i - 1); } } if (Blk[6] && nsol >= 1) dprtsv( a, mda, n, nsol, names,names_s, 2, unit, width ); return; #undef NAMES #undef A } /* end of function */