/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:40 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "daccum.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	MARK	2
#define	ZERO	0.0e0
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ daccum(
double *a,
long lda,
long n,
double *b,
long ldb,
long nb,
long *ir1,
long nrows,
long *ncount)
{
#define A(I_,J_)	(*(a+(I_)*(lda)+(J_)))
#define B(I_,J_)	(*(b+(I_)*(ldb)+(J_)))
	long int i, ip, irow, j, l1, last, m;
	double cfac, sfac, temp, uparam;
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 1996-03-30 DACCUM Krogh   Added external statement.
	 *>> 1994-11-11 DACCUM Krogh   Declared all vars.
	 *>> 1994-10-20 DACCUM Krogh  Changes to use M77CON
	 *>> 1987-11-24 DACCUM Lawson  Initial code.
	 *--D replaces "?": ?ACCUM, ?ROTG, ?HTCC, ?NRM2
	 *     Sequential accumulation of rows of data for a linear least
	 *     squares problem.  Using Givens orthogonal transformations when
	 *     NROWS is small and Householder orthogonal transformations when
	 *     NROWS is larger.
	 *
	 *        On first call user must set IR1 = 1.  On each return this subr
	 *     will update IR1 to min(IR1+NROWS, N+2).  The user should not alter
	 *     IR1 while processing data for the same case.
	 *
	 *        On all calls, including the first, the user sets NROWS to
	 *     indicate the number of rows of new data being provided.  The user
	 *     must put the new data in rows IR1 through IR1 + NROWS - 1 of the
	 *     arrays A() and B().  On return, with IR1 updated to
	 *     min(IR1+NROWS, N+2), the first IR1-1 rows of A() will contain
	 *     transformed data on and to the right of the diagonal, and zeros
	 *     to the left of the diagonal.  The first IR1-1 rows of B() will
	 *     also contain transformed data.
	 *     ------------------------------------------------------------------
	 *     Reference: C. L. Lawson & R. J. Hanson,
	 *                Solving Least Squares Problems, Prentice-Hall, 1974.
	 *     Original code by R. J. Hanson, JPL, Sept 11, 1968.
	 *     Adapted to Fortran 77 for the JPL MATH77 library
	 *     by Lawson, 6/2/87.
	 *     ------------------------------------------------------------------
	 *                      Subroutine arguments
	 *
	 *     A(,)  [inout]
	 *     LDA   [in]
	 *     N     [in]
	 *     B(,)  [inout]
	 *     LDB   [in]
	 *     NB    [in]
	 *     IR1   [inout]
	 *     NROWS  [in]
	 *     NCOUNT  [out]
	 *     ------------------------------------------------------------------
	 *     Subprograms called directly: DHTCC, DNRM2, DROTG, IERM1, IERV1
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------
	 * */
	if (nrows <= 0)
		return;
 
	m = *ir1 + nrows - 1;
	if (m > lda || m > ldb)
	{
		ierm1( "DACCUM", 1, 0, "Require LDA .ge. M and LDB .ge M where M = IR1+NROWS-1"
		 , "IR1", *ir1, ',' );
		ierv1( "NROWS", nrows, ',' );
		ierv1( "LDA", lda, ',' );
		ierv1( "LDB", ldb, '.' );
		return;
	}
	if (m <= n)
	{
		last = m - 1;
	}
	else
	{
		last = n;
	}
 
	if (nrows <= MARK)
	{
		/*      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
		 *           Use Givens transformation for better efficiency when the
		 *           number of rows being introduced is .le. MARK.
		 * */
		for (ip = 1; ip <= last; ip++)
		{
			l1 = max( *ir1, ip + 1 );
			for (irow = l1; irow <= m; irow++)
			{
				drotg( &A(ip - 1,ip - 1), &A(ip - 1,irow - 1), &cfac,
				 &sfac );
				A(ip - 1,irow - 1) = ZERO;
				for (j = ip + 1; j <= n; j++)
				{
					temp = A(j - 1,ip - 1);
					A(j - 1,ip - 1) = cfac*temp + sfac*A(j - 1,irow - 1);
					A(j - 1,irow - 1) = -sfac*temp + cfac*A(j - 1,irow - 1);
				}
				for (j = 1; j <= nb; j++)
				{
					temp = B(j - 1,ip - 1);
					B(j - 1,ip - 1) = cfac*temp + sfac*B(j - 1,irow - 1);
					B(j - 1,irow - 1) = -sfac*temp + cfac*B(j - 1,irow - 1);
				}
			}
		}
		/*        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
	}
	else
	{
		/*        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
		 *                 Use Householder transformations
		 *                 when the number of rows being introduced exceeds MARK.
		 * */
		for (ip = 1; ip <= last; ip++)
		{
			l1 = max( *ir1, ip + 1 );
			dhtcc( 1, ip, l1, m, &A(ip - 1,0), &uparam, &A(min( ip + 1, n ) - 1,0),
			 lda, n - ip );
			dhtcc( 2, ip, l1, m, &A(ip - 1,0), &uparam, &B(0,0), ldb,
			 nb );
 
			/*                      Clear elements just made implicitly zero.
			 * */
			for (i = l1; i <= min( m, n + 1 ); i++)
			{
				A(ip - 1,i - 1) = ZERO;
			}
		}
		/*        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
	}
 
	if (*ir1 == 1)
	{
		*ncount = nrows;
	}
	else
	{
		*ncount += nrows;
	}
	*ir1 = min( *ir1 + nrows, n + 2 );
	if (m <= n + 1)
		return;
 
	/*     Pack lengths of B() column vectors below row N
	 *     to single locations each.
	 * */
	for (j = 1; j <= nb; j++)
	{
		B(j - 1,n) = dnrm2( m - n, &B(j - 1,n), 1 );
	}
	return;
#undef	B
#undef	A
} /* end of function */