/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:45 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sbacc.h" #include /* PARAMETER translations */ #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ sbacc( float *g, long ldg, long nb, long *ir, long mt, long jt, long *jtprev, long *ierr2) { #define G(I_,J_) (*(g+(I_)*(ldg)+(J_))) long int i, ie, ierold, ig, ig1, ig2, istep, j, jg, k, kh, kh1, l, lp1, mh, nbp1; float uparam; /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1994-10-20 SBACC Krogh Changes to use M77CON *>> 1992-06-16 SBACC CLL *>> 1989-10-20 CLL *>> 1987-11-24 SBACC Lawson Initial code. * Sequential accumulation of data for a banded least-squares * problem. For solution use _BSOL. * ------------------------------------------------------------------ * Subroutine arguments * * G(,) [inout] Array in which user provides data and this subr * builds a packed representation of the orthogonally * transformed equivalent problem. * LDG [in] Leading dimensioning parameter for G(,). * NB [in] Bandwidth of data matrix. Set on first call for a * new problem and not changed thereafter. * IR [inout] User must set IR = 1 on 1st call for a new problem. * User must not thereafter alter IR. This subr will update * IR as IR := JT + min(NB+1, MT+max(IR-JT,0)). * IR identifies the row of G(,) in which user's new block * of data begins. * MT [in] Number of new rows of data being provided on the * current entry. * JT [in] Column index in the user's problem of the data being * provided in column 1 of G(,). * JTPREV [inout] Should not be set or altered by the user. * This subr will update JTPREV as JTPREV := JT. * JTPREV identifies a row in G(,) at which the mode of * packed storage changes. Data in row JTPREV and below * are subject to possible change as future data are * introduced, whereas data above row JTPREV are not. * IERR2 [inout] Error status indicator. Must be set to zero on * initial call to this subr for a new problem. Values on * return have meaning as follows: * 0 means no errors detected. * 1 means MT exceeds LDG - IR + 1 * 2 means JT .lt. JTPREV * 3 means JT .gt. min(JTPREV + NB, IR) which causes the * problem to be structurally singular. The transformations * by this subr can be completed in this case but the * resulting transformed problem cannot be solved simply * by use of _BSOL. * 4 means the condition of Error 3 above applys but computa- * tion must be abandonded due to storage limitations * indicated by MT > LDG - JT + 1. * 5 This subr will STOP because the input value of IERR2 was * nonzero. * ------------------------------------------------------------------ * C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 * Reference: 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974 * Comments in this code refer to Algorithm Step numbers on * pp 213-217 of the book. * ------------------------------------------------------------------ * 1984 July 11, C. L. Lawson, JPL. Adapted code from book * to Fortran 77 for JPL MATH 77 library. * Prefixing subprogram names with S or D for s.p. or d.p. versions. * Using generic names for any intrinsic functions. * July 1987. Using _HTCC in place of _H12. * 1989-10-20 CLL Moved integer declaration earlier to avoid warning * msg from Cray compiler. * 1992-06-16 CLL Using "min" fcn in arg list of call to _HTCC * to avoid "index out of range" when using a bounds checker, even * though no reference would be made to the "out of range" location. * ------------------------------------------------------------------ *--S replaces "?": ?BACC, ?HTCC * Both versions use IERM1, IERV1 * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * * ALG. STEPS 1-4 ARE PERFORMED EXTERNAL TO THIS SUBROUTINE. * */ if (mt <= 0) return; if (*ir == 1) { /* Initialize for new problem. */ *jtprev = 1; *ierr2 = 0; } else if (*ierr2 != 0) { ierold = *ierr2; *ierr2 = 5; ierm1( "SBACC", *ierr2, 0, "Value of IERR2 nonzero on entry" , "IERR2 on entry", ierold, '.' ); exit(0); } nbp1 = nb + 1; if (mt > ldg - *ir + 1) { *ierr2 = 1; ierm1( "SBACC", *ierr2, 0, "MT must not exceed LDG - IR + 1" , "MT", mt, ',' ); ierv1( "LDG", ldg, ',' ); ierv1( " IR", *ir, '.' ); return; } /* ALG. STEP 5 */ if (jt == *jtprev) goto L_70; /* ALG. STEP 6 */ if (jt < *jtprev) { *ierr2 = 2; ierm1( "SBACC", *ierr2, 0, "Require JT .ge. JTPREV", "JT", jt, ',' ); ierv1( "JTPREV", *jtprev, '.' ); return; } /* ALG. STEP 7 */ if (jt > min( *jtprev + nb, *ir )) { *ierr2 = 3; ierm1( "SBACC", *ierr2, 0, "With JT .gt. min(JTPREV+NB,IR) problem will be singular." , "JT", jt, ',' ); ierv1( "JTPREV", *jtprev, ',' ); ierv1( "NB", nb, ',' ); ierv1( "IR", *ir, '.' ); /* Note: At this point we know the problem will be structurally * singular and thus will not be solvable simply by use of _BSOL. * This does not preclude, however, completion of the transforma- * tion to the band triangular form and so we continue. The * calling program will be informed of this condition by the set- * ting of of IERR2 = 3, and thus can quit or continue as desired. * If this is not due to a usage mistake * then the user should probably introduce more data, modify his * or her mathematical model, or use a stabilizing method such as * is discussed on pp. 218-219 of the L & H book. * */ if (mt > ldg - jt + 1) { *ierr2 = 4; ierm1( "SBACC", *ierr2, 0, "MT must not exceed LDG - JT + 1" , "MT", mt, ',' ); ierv1( "LDG", ldg, ',' ); ierv1( "JT", jt, '.' ); return; } /* ALG. STEPS 8-9 */ for (i = 1; i <= mt; i++) { ig1 = jt + mt - i; ig2 = *ir + mt - i; for (j = 1; j <= nbp1; j++) { G(j - 1,ig1 - 1) = G(j - 1,ig2 - 1); } } /* ALG. STEP 10 */ ie = jt - *ir; for (i = 1; i <= ie; i++) { ig = *ir + i - 1; for (j = 1; j <= nbp1; j++) { G(j - 1,ig - 1) = ZERO; } } /* ALG. STEP 11 */ *ir = jt; } /* Here we must shift some rows left by variable amounts. * min(NB-1,IR-JTPREV-1) is the number of rows to be shifted. * K = min(L,JT-JTPREV) is the extent of the left shift * for row JTPREV+L. * * ALG. STEP 12-13 */ for (l = 1; l <= min( nb - 1, *ir - *jtprev - 1 ); l++) { /* ALG. STEP 14 */ k = min( l, jt - *jtprev ); /* ALG. STEP 15 */ lp1 = l + 1; ig = *jtprev + l; for (i = lp1; i <= nb; i++) { G(i - k - 1,ig - 1) = G(i - 1,ig - 1); } /* ALG. STEP 16 */ for (jg = nbp1 - k; jg <= nb; jg++) { G(jg - 1,ig - 1) = ZERO; } } /* ALG. STEP 17 */ *jtprev = jt; /* ALG. STEPS 18-19 */ L_70: ; /* MH = No. of rows in the block to which the orthogonal * triangularization will be applied. * KH1 = Number of pivot columns to be used in the orghogonal * transformations. * */ mh = *ir + mt - *jtprev; kh1 = min( nbp1, mh - 1 ); istep = *ir - *jtprev + 1; /* ALG. STEP 20 */ for (i = 1; i <= kh1; i++) { shtcc( 1, i, max( i + 1, istep ), mh, &G(i - 1,*jtprev - 1), &uparam, &G(min( i + 1, nbp1 ) - 1,*jtprev - 1), ldg, nbp1 - i ); } /* ALG. STEP 21 */ kh = min( nbp1, mh ); *ir = *jtprev + kh; /* ALG. STEP 22 */ if (kh >= nbp1) { /* ALG. STEP 23 */ for (i = 1; i <= nb; i++) { G(i - 1,*ir - 2) = ZERO; } /* ALG. STEP 24 */ } /* ALG. STEP 25 */ return; #undef G } /* end of function */