/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:17 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dqrbd.h" #include #include /* PARAMETER translations */ #define ONE 1.0e0 #define TEN 10.0e0 #define TWO 2.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dqrbd( long *ipass, double q[], double e[], long n, double *v, long mdv, long nrv, double *c, long mdc, long ncc) { #define V(I_,J_) (*(v+(I_)*(mdv)+(J_))) #define C(I_,J_) (*(c+(I_)*(mdc)+(J_))) LOGICAL32 fail, havers, wntv; long int _l0, i, ii, j, k, l, lp1, n10, nqrs; double big, cs, dnorm, eps, f, g, h, small, sn, t, x, y, z; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const E = &e[0] - 1; double *const Q = &q[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-03-30 DQRBD Krogh Added external statement. *>> 1994-10-20 DQRBD Krogh Changes to use M77CON *>> 1992-03-13 DQRBD FTK Removed implicit statements. *>> 1987-11-24 DQRBD Lawson Initial code. *--D replaces "?": ?QRBD, ?ROT, ?ROTG, ?SWAP * * Computes the singular value decomposition of an N order * upper bidiagonal matrix, B. This decomposition is of the form * * (1) B = U * S * (V**t) * * where U and V are each N x N orthogonal and S is N x N * diagonal, with nonnegative diagonal terms in nonincreasing * order. Note that these matrices also satisfy * * S = (U**t) * B * V * * The user may optionally provide an NRV x N matrix V1 in * the array V(,) and an N x NCC matrix C in the array C(,). * This subr will replace these matrices respectively by the * NRV x N matrix * V2 = V1 * V * and the N x NCC matrix * C2 = (U**t) * C1 * * On entry the bidiagonal matrix, B, is to be given in the * arrays Q() and E() as indicated by the diagram: * * (Q1,E2,0... ) * ( Q2,E3,0... ) * B = ( . ) * ( . 0) * ( 0 .EN) * ( QN) * * On return, the singular values, i.e. the diagonal terms of * S, will be stored in Q(1) through Q(N). * * THIS CODE IS BASED ON THE PAPER AND 'ALGOL' CODE.. * REF.. * 1. REINSCH,C.H. AND GOLUB,G.H. 'SINGULAR VALUE DECOMPOSITION * AND LEAST SQUARES SOLUTIONS' (NUMER. MATH.), VOL. 14,(1970). * * ------------------------------------------------------------------ * Subroutine Arguments * * IPASS (Out) On return, IPASS = 1 means computation was * successful or N .le. 0. * IPASS = 2 means computation not successful after 10*N * iterations. The algorithm usually succeeds in about * 2*N iterations. * * Q() (In/Out) On entry must contain the diagonal terms of * the bidiagonal matrix B. On return contains the N * singular values of B. These will be nonnegative and * in nonincreasing order. * * E() (In/Work) On entry must contain the superdiagonal terms * of B. E(1) is not used. For i = 2, ..., N, E(i) must * contain B(i-1,i). The contents of this array will be * modified by this subr. * * N (In) The order of the bidiagonal matrix B and the * number of singular values to be produced. * * V(,) (In/Out) If NRV .gt. 0, on entry this array contains an * NRV x N matrix, V1. On return V(,) will contain the * NRV x N matrix V2 = V1 * V where V is defined by Eq.(1) * above. If NRV .eq. 0 the array V(,) will not be * referenced. * * MDV (In) First dimensioning parameter for the array V(,). * Require MDV .ge. NRV. * * NRV (In) No. of rows in matrix V1 in array V(,). * Require NRV .ge. 0. * * C(,) (In/Out) If NCC .gt. 0, on entry this array contains an * N x NCC matrix, C1. On return C(,) will contain the * N x NCC matrix C2 = (U**t) * C1 where U is defined by * Eq.(1) above. If NCC .eq. 0 the array C(,) will not be * referenced. * * MDC (In) First dimensioning parameter for the array C(,). * Require MDC .ge. N. * * NCC (In) No. of columns in matrix C1 in array C(,). * Require NCC .ge. 0. * ------------------------------------------------------------------ * Subprograms referenced: ERMSG, D1MACH, DROT, DROTG, DSWAP. * ------------------------------------------------------------------ * 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, 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. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ eps = DBL_EPSILON; big = TEN/sqrt( eps ); *ipass = 1; if (n <= 0) return; n10 = 10*n; wntv = nrv > 0; havers = ncc > 0; fail = FALSE; nqrs = 0; E[1] = ZERO; dnorm = ZERO; for (j = 1; j <= n; j++) { dnorm = fmax( fabs( Q[j] ) + fabs( E[j] ), dnorm ); } small = dnorm*eps; /* ------------------------------------------------------------------ * Begin main loop. */ for (k = n; k >= 1; k--) { /* TEST FOR SPLITTING OR RANK DEFICIENCIES.. * FIRST MAKE TEST FOR LAST DIAGONAL TERM, Q(K), BEING SMALL. */ L_20: ; if (k == 1) goto L_50; if (fabs( Q[k] ) > small) goto L_50; /* SINCE Q(K) IS SMALL WE WILL MAKE A SPECIAL PASS TO * TRANSFORM E(K) TO ZERO. * */ cs = ZERO; sn = -ONE; for (ii = 2; ii <= k; ii++) { i = k + 1 - ii; f = -sn*E[i + 1]; E[i + 1] *= cs; drotg( &Q[i], &f, &cs, &sn ); /* TRANSFORMATION CONSTRUCTED TO ZERO POSITION (I,K). * * ACCUMULATE RT. TRANSFORMATIONS IN V. */ if (wntv) drot( nrv, &V(i - 1,0), 1, &V(k - 1,0), 1, cs, sn ); } /* THE MATRIX IS NOW BIDIAGONAL, AND OF LOWER ORDER * SINCE E(K) .EQ. ZERO.. * */ L_50: for (l = k; l >= 1; l--) { if (fabs( E[l] ) <= small) goto L_100; if (fabs( Q[l - 1] ) <= small) goto L_70; } /* THE ABOVE LOOP CAN'T COMPLETE SINCE E(1) = ZERO. * */ goto L_100; /* CANCELLATION OF E(L), L.GT.1. * */ L_70: cs = ZERO; sn = -ONE; for (i = l; i <= k; i++) { f = -sn*E[i]; E[i] *= cs; if (fabs( f ) <= small) goto L_100; drotg( &Q[i], &f, &cs, &sn ); if (havers) drot( ncc, &C(0,i - 1), mdc, &C(0,l - 2), mdc, cs, sn ); } /* TEST FOR CONVERGENCE.. * */ L_100: z = Q[k]; if (l == k) goto L_170; /* SHIFT FROM BOTTOM 2 BY 2 MINOR OF B**(T)*B. * */ x = Q[l]; y = Q[k - 1]; g = E[k - 1]; h = E[k]; f = ((y - z)*(y + z) + (g - h)*(g + h))/(TWO*h*y); if (fabs( f ) < big) { g = sqrt( ONE + SQ(f) ); } else { g = fabs( f ); } if (f >= ZERO) { t = f + g; } else { t = f - g; } f = ((x - z)*(x + z) + h*(y/t - h))/x; /* NEXT QR SWEEP.. */ cs = ONE; sn = ONE; lp1 = l + 1; for (i = lp1; i <= k; i++) { g = E[i]; y = Q[i]; h = sn*g; g *= cs; drotg( &f, &h, &cs, &sn ); E[i - 1] = f; f = x*cs + g*sn; g = -x*sn + g*cs; h = y*sn; y *= cs; /* ACCUMULATE ROTATIONS (FROM THE RIGHT) IN 'V' * */ if (wntv) drot( nrv, &V(i - 2,0), 1, &V(i - 1,0), 1, cs, sn ); drotg( &f, &h, &cs, &sn ); Q[i - 1] = f; f = cs*g + sn*y; x = -sn*g + cs*y; /* APPLY ROTATIONS FROM THE LEFT TO RIGHT HAND SIDES IN 'C' * */ if (havers) drot( ncc, &C(0,i - 2), mdc, &C(0,i - 1), mdc, cs, sn ); } E[l] = ZERO; E[k] = f; Q[k] = x; nqrs += 1; if (nqrs <= n10) { /* Return for further iteration on Q(K). */ goto L_20; } else { /* Iteration count limit exceeded for Q(K). */ fail = TRUE; } /* ------------------------------------------------------------------ * Accepting Q(K). Either convergence has been reached or * the iteration count limit has been exceeded. * Adjust sign of Q(K) to be nonnegative. * */ L_170: if (z < ZERO) { Q[k] = -z; if (wntv) { for (j = 1; j <= nrv; j++) { V(k - 1,j - 1) = -V(k - 1,j - 1); } } } } /* End of main loop. * ------------------------------------------------------------------ */ if (n == 1) return; /* Test sing vals for being in order. If not in order then * sort them. * */ for (i = 2; i <= n; i++) { if (Q[i] > Q[i - 1]) goto L_220; } /* Sing vals are sorted. */ goto L_290; /* Sort the sing vals. */ L_220: for (i = 2; i <= n; i++) { t = Q[i - 1]; k = i - 1; for (j = i; j <= n; j++) { if (t < Q[j]) { t = Q[j]; k = j; } } if (k != i - 1) { Q[k] = Q[i - 1]; Q[i - 1] = t; if (havers) dswap( ncc, &C(0,i - 2), mdc, &C(0,k - 1), mdc ); if (wntv) dswap( nrv, &V(i - 2,0), 1, &V(k - 1,0), 1 ); } } /* Sing vals are sorted. * * ------------------------------------------------------------------ */ L_290: ; if (fail) { /* Error message. */ *ipass = 2; ermsg( "DQRBD", 1, 0, "Convergence failure in computing singular values." , '.' ); } return; #undef V #undef C } /* end of function */