/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:55 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sq7rfh.h" #include #include /* PARAMETER translations */ #define ONE 1.0e0 #define TEN 1.e1 #define WTOL 0.75e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ sq7rfh( long *ierr, long ipivot[], long n, long nn, long nopivk, long p, float *q, float r[], long rlen, float w[]) { #define Q(I_,J_) (*(q+(I_)*(nn)+(J_))) long int i, ii, j, k, kk, km1, kp1, nk1; float ak, big, qkk, s, singtl, t, t1, wk; static float bigrt = 0.0e0; static float meps10 = 0.0e0; static float tiny = 0.e0; static float tinyrt = 0.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Ipivot = &ipivot[0] - 1; float *const R = &r[0] - 1; float *const W = &w[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. * File: SQ7RFH.for Four subrs used by the * David Gay & Linda Kaufman nonlinear LS package. * SQ7RFH, SV7PRM, SS7CPR, SV7SWP * Needed by versions that handle the Separable problem. * The first 3 are called by DRNSG & DRNSGB. * SV7SWP is called by SQ7RFH. * *>> 1998-10-29 SQ7RFH Krogh Moved external statement up for mangle. *>> 1996-05-15 SQ7RFH Krogh Use M77CON for conversion to C. *>> 1995-11-15 SQ7RFH Krogh Moved formats up for C conversion. *>> 1994-11-02 SQ7RFH Krogh Changes to use M77CON *>> 1992-04-27 SQ7RFH CLL Removed unreferenced stmt labels. *>> 1992-04-13 CLL Change from Hollerith to '...' syntax in formats. *>> 1990-03-30 CLL @ JPL *>> 1990-02-21 CLL @ JPL *--S replaces "?": ?Q7RFH,?D7TPR,?R7MDC,?V2AXY,?V2NRM,?V7SCL,?V7SCP *--& ?V7SWP,?S7CPR,?V7PRM * * *** COMPUTE QR FACTORIZATION VIA HOUSEHOLDER TRANSFORMATIONS * *** WITH COLUMN PIVOTING *** * * *** PARAMETER DECLARATIONS *** * */ /* DIMENSION R(P*(P+1)/2) * * --------------------------- DESCRIPTION ---------------------------- * * THIS ROUTINE COMPUTES A QR FACTORIZATION (VIA HOUSEHOLDER TRANS- * FORMATIONS) OF THE MATRIX A THAT ON INPUT IS STORED IN Q. * IF NOPIVK ALLOWS IT, THIS ROUTINE DOES COLUMN PIVOTING -- IF * K .GT. NOPIVK, THEN ORIGINAL COLUMN K IS ELIGIBLE FOR PIVOTING. * THE Q AND R RETURNED ARE SUCH THAT COLUMN I OF Q*R EQUALS * COLUMN IPIVOT(I) OF THE ORIGINAL MATRIX A. THE UPPER TRIANGULAR * MATRIX R IS STORED COMPACTLY BY COLUMNS, I.E., THE OUTPUT VECTOR R * CONTAINS R(1,1), R(1,2), R(2,2), R(1,3), R(2,3), ..., R(P,P) (IN * THAT ORDER). IF ALL GOES WELL, THEN THIS ROUTINE SETS IERR = 0. * BUT IF (PERMUTED) COLUMN K OF A IS LINEARLY DEPENDENT ON * (PERMUTED) COLUMNS 1,2,...,K-1, THEN IERR IS SET TO K AND THE R * MATRIX RETURNED HAS R(I,J) = 0 FOR I .GE. K AND J .GE. K. * THE ORIGINAL MATRIX A IS AN N BY P MATRIX. NN IS THE LEAD * DIMENSION OF THE ARRAY Q AND MUST SATISFY NN .GE. N. NO * PARAMETER CHECKING IS DONE. * PIVOTING IS DONE AS THOUGH ALL COLUMNS OF Q WERE FIRST * SCALED TO HAVE THE SAME NORM. IF COLUMN K IS ELIGIBLE FOR * PIVOTING AND ITS (SCALED) NORM**2 LOSS IS MORE THAN THE * MINIMUM SUCH LOSS (OVER COLUMNS K THRU P), THEN COLUMN K IS * SWAPPED WITH THE COLUMN OF LEAST NORM**2 LOSS. * * CODED BY DAVID M. GAY (FALL 1979, SPRING 1984). * * ------------------------- LOCAL VARIABLES -------------------------- * */ /*+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ * */ *ierr = 0; if (meps10 > ZERO) goto L_10; bigrt = sr7mdc( 5 ); meps10 = TEN*sr7mdc( 3 ); tinyrt = sr7mdc( 2 ); tiny = sr7mdc( 1 ); big = sr7mdc( 6 ); if (tiny*big < ONE) tiny = ONE/big; L_10: singtl = (float)( max( n, p ) )*meps10; /* *** INITIALIZE W, IPIVOT, AND DIAG(R) *** * */ j = 0; for (i = 1; i <= p; i++) { Ipivot[i] = i; t = sv2nrm( n, &Q(i - 1,0) ); if (t > ZERO) goto L_20; W[i] = ONE; goto L_30; L_20: W[i] = ZERO; L_30: j += i; R[j] = t; } /* *** MAIN LOOP *** * */ kk = 0; nk1 = n + 1; for (k = 1; k <= p; k++) { if (nk1 <= 1) goto L_999; nk1 -= 1; kk += k; kp1 = k + 1; if (k <= nopivk) goto L_60; if (k >= p) goto L_60; /* *** FIND COLUMN WITH MINIMUM WEIGHT LOSS *** * */ t = W[k]; if (t <= ZERO) goto L_60; j = k; for (i = kp1; i <= p; i++) { if (W[i] >= t) goto L_50; t = W[i]; j = i; L_50: ; } if (j == k) goto L_60; /* *** INTERCHANGE COLUMNS K AND J *** * */ i = Ipivot[k]; Ipivot[k] = Ipivot[j]; Ipivot[j] = i; W[j] = W[k]; W[k] = t; i = j*(j + 1)/2; t1 = R[i]; R[i] = R[kk]; R[kk] = t1; sv7swp( n, &Q(k - 1,0), &Q(j - 1,0) ); if (k <= 1) goto L_60; i += -j + 1; j = kk - k + 1; sv7swp( k - 1, &R[i], &R[j] ); /* *** COLUMN K OF Q SHOULD BE NEARLY ORTHOGONAL TO THE PREVIOUS * *** COLUMNS. NORMALIZE IT, TEST FOR SINGULARITY, AND DECIDE * *** WHETHER TO REORTHOGONALIZE IT. * */ L_60: ak = R[kk]; if (ak <= ZERO) goto L_140; wk = W[k]; /* *** SET T TO THE NORM OF (Q(K,K),...,Q(N,K)) * *** AND CHECK FOR SINGULARITY. * */ if (wk < WTOL) goto L_70; t = sv2nrm( nk1, &Q(k - 1,k - 1) ); if (t/ak <= singtl) goto L_140; goto L_80; L_70: t = sqrtf( ONE - wk ); if (t <= singtl) goto L_140; t *= ak; /* *** DETERMINE HOUSEHOLDER TRANSFORMATION *** * */ L_80: qkk = Q(k - 1,k - 1); if (t <= tinyrt) goto L_90; if (t >= bigrt) goto L_90; if (qkk < ZERO) t = -t; qkk += t; s = sqrtf( t*qkk ); goto L_110; L_90: s = sqrtf( t ); if (qkk < ZERO) goto L_100; qkk += t; s *= sqrtf( qkk ); goto L_110; L_100: t = -t; qkk += t; s *= sqrtf( -qkk ); L_110: Q(k - 1,k - 1) = qkk; /* *** SCALE (Q(K,K),...,Q(N,K)) TO HAVE NORM SQRT(2) *** * */ if (s <= tiny) goto L_140; sv7scl( nk1, &Q(k - 1,k - 1), ONE/s, &Q(k - 1,k - 1) ); R[kk] = -t; /* *** COMPUTE R(K,I) FOR I = K+1,...,P AND UPDATE Q *** * */ if (k >= p) goto L_999; j = kk + k; ii = kk; for (i = kp1; i <= p; i++) { ii += i; sv2axy( nk1, &Q(i - 1,k - 1), -sd7tpr( nk1, &Q(k - 1,k - 1), &Q(i - 1,k - 1) ), &Q(k - 1,k - 1), &Q(i - 1,k - 1) ); t = Q(i - 1,k - 1); R[j] = t; j += i; t1 = R[ii]; if (t1 > ZERO) W[i] += powif(t/t1,2); } } /* *** SINGULAR Q *** * */ L_140: *ierr = k; km1 = k - 1; j = kk; for (i = k; i <= p; i++) { sv7scp( i - km1, &R[j], ZERO ); j += i; } L_999: return; /* *** LAST CARD OF SQ7RFH FOLLOWS *** */ #undef Q } /* end of function */ /* PARAMETER translations */ #define PRUNIT 21 #define SOLPRT 22 /* end of PARAMETER translations */ void /*FUNCTION*/ ss7cpr( float c[], long iv[], long l, long liv) { long int pu; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const C = &c[0] - 1; long *const Iv = &iv[0] - 1; /* end of OFFSET VECTORS */ /* *** PRINT C FOR DNSG (ETC.) *** * * ------------------------------------------------------------------ */ long int i; /* INTEGER I * * */ /* *** BODY *** * */ if (Iv[1] > 11) goto L_999; if (Iv[SOLPRT] == 0) goto L_999; pu = Iv[PRUNIT]; if (pu == 0) goto L_999; /*++ CODE for .C. is active */ if( l > 0 ) { printf("\n LINEAR PARAMETERS...\n"); for( i = 1; i <= l; i++ ) printf("\n%5ld%16.6g", i, C[i]); printf("\n");} L_999: return; /*++ CODE for ~.C. is inactive * IF (L .GT. 0) WRITE(PU,10) (I, C(I), I = 1, L) *++ END * * *** LAST LINE OF SS7CPR FOLLOWS *** */ } /* end of function */ void /*FUNCTION*/ sv7prm( long n, long ip[], float x[]) { long int i, j, k; float s, t; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Ip = &ip[0] - 1; float *const X = &x[0] - 1; /* end of OFFSET VECTORS */ /* PERMUTE X SO THAT X.OUTPUT(IP(I)) = X.INPUT(I). * IP IS UNCHANGED ON OUTPUT. * * ------------------------------------------------------------------ */ for (i = 1; i <= n; i++) { j = Ip[i]; if (j == i) goto L_30; if (j > 0) goto L_10; Ip[i] = -j; goto L_30; L_10: t = X[i]; L_20: s = X[j]; X[j] = t; t = s; k = j; j = Ip[k]; Ip[k] = -j; if (j > i) goto L_20; X[j] = t; L_30: ; } /* *** LAST LINE OF SV7PRM FOLLOWS *** */ return; } /* end of function */ void /*FUNCTION*/ sv7swp( long n, float x[], float y[]) { long int i; float t; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const X = &x[0] - 1; float *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* *** INTERCHANGE N-VECTORS X AND Y. *** * * ------------------------------------------------------------------ */ for (i = 1; i <= n; i++) { t = X[i]; X[i] = Y[i]; Y[i] = t; } /* *** LAST CARD OF SV7SWP FOLLOWS *** */ return; } /* end of function */