/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:54 */
/*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 "sdasco.h"
#include <stdlib.h>
void /*FUNCTION*/ sdasco(
float *a,
long lda,
long neq,
float wt[],
float delta[])
{
#define A(I_,J_)	(*(a+(I_)*(lda)+(J_)))
	long int i, imax, j, nc, nrank;
	float c, s, tn, z;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const Delta = &delta[0] - 1;
	float *const Wt = &wt[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 2006, Math a la Carte, Inc.
	 *>> 2006-05-20 sdasco Hanson Fix bug for computing norm
	 *>> 2006-05-18 sdasco Hanson return residual norm value in delta(neq+1)
	 *>> 2006-04-26 sdasco Hanson installed pivoting and test for def rank
	 *>> 2003-03-06 sdasco Hanson accounted for use of reciprocal weights
	 *>> 2002-03-19 sdasco Krogh  Added sqrt on the weights.
	 *>> 2001-11-23 sdasco Krogh  Changed many names per library conventions.
	 *>> 2001-11-04 sdasco Krogh  Fixes for F77 and conversion to single
	 *>> 2001-11-01 sdasco Hanson Provide code to Math a la Carte.
	 *--S replaces "?":?dasco, ?nrm2, ?dot, ?scal, ?rotg, ?rot, ?copy, ?swap,
	 *--& ?axpy, i?amax
	 *      implicit none
	 *     Compute a single Newton step:
	 *     Solve an underdetermined system.
	 *     The number of constraint rows is NC=LDD-NEQ.
	 *
	 *     On input the first neq entries of delta contain the reciprocals of
	 *     the weights.  The next entries contain the constraint residuals.
	 *     On the other hand, the matrix, a, contains only the partials for
	 *     the constraint rows.
	 * */
	nc = lda - neq;
	if (nc <= 0)
		return;
	/*     Use weights derived from error norm. */
	for (j = 1; j <= neq; j++)
	{
		sscal( nc, sqrtf( Wt[j] ), &A(j - 1,0), 1 );
	}
 
	/*     Triangularize constraint matrix using right-handed
	 *     plane rotations.  These are stored in the entry
	 *     eliminated and then reconstructed. */
	for (i = 1; i <= min( neq, nc - 1 ); i++)
	{
		/*     Do row pivoting to maximize diagonal term (in L2 norm)
		 *     so that smallest pivots are at the SE part of the lower triangle. */
		for (j = i; j <= nc; j++)
		{
			Delta[j] = sdot( neq - i + 1, &A(i - 1,j - 1), lda, &A(i - 1,j - 1),
			 lda );
		}
		imax = isamax( nc - i + 1, &Delta[i], 1 );
		/*     If found a row with larger L2 norm, interchange this row
		 *     and the right-hand side. */
		if (imax > 1)
		{
			sswap( neq, &A(0,i - 1), lda, &A(0,i + imax - 2), lda );
			sswap( 1, &Delta[neq + i], 1, &Delta[neq + i + imax - 1],
			 1 );
		}
		/*     Eliminate the non-zero entries and store the rotation
		 *     in the place eliminated.  These are used to compute
		 *     the least-distance move back to the constraints.
		 *     Note that this uses an often overlooked feature of ?rotg. */
		for (j = i + 1; j <= neq; j++)
		{
			srotg( &A(i - 1,i - 1), &A(j - 1,i - 1), &c, &s );
			srot( nc - i, &A(i - 1,i), 1, &A(j - 1,i), 1, c, s );
		}
	}
	/*     This is a special elimination step for the last row. */
	for (j = nc + 1; j <= neq; j++)
	{
		srotg( &A(nc - 1,nc - 1), &A(j - 1,nc - 1), &c, &s );
	}
	/*     System is now lower triangular.  Forward solve step.
	 *     The right-hand side for this system is in DELTA(NEQ+I),
	 *     I=1,NC. */
	Delta[1] = 0.e0;
	scopy( neq, delta, 0, delta, 1 );
	nrank = min( nc, neq );
	for (i = 1; i <= nc; i++)
	{
		/*     Due to the weighting used, a small number is .le. 1 for
		 *     a diagonal term. */
		if (fabsf( A(i - 1,i - 1) ) <= 1.e0)
		{
			nrank = i - 1;
			goto L_55;
		}
		Delta[i] = Delta[neq + i];
	}
L_55:
	;
	/* Solve for the update. */
	for (i = 1; i <= (nrank - 1); i++)
	{
		z = Delta[i]/A(i - 1,i - 1);
		saxpy( nrank - i, -z, &A(i - 1,i), 1, &Delta[i + 1], 1 );
		Delta[i] = z;
	}
	Delta[nrank] /= A(nrank - 1,nrank - 1);
	/*     Compute residuals on those equations not solved.
	 *     They should be consistent but may not due to
	 *     modelling errors or large numerical errors. */
	for (i = nrank + 1; i <= nc; i++)
	{
		Delta[neq + i] -= sdot( nrank, &A(0,i - 1), lda, delta, 1 );
	}
	/*     Now apply the Givens transformations in reverse order.
	 *     They are stored in one entry form so the values are
	 *     reconstructed as in the 1979 Level 1 BLAS paper or
	 *     LINPACK, page A.5. */
	for (i = min( neq, nc ); i >= 1; i--)
	{
		for (j = neq; j >= (i + 1); j--)
		{
			z = A(j - 1,i - 1);
			if (fabsf( z ) < 1.e0)
			{
				s = z;
				c = sqrtf( (1.e0 - s)*(1.e0 + s) );
			}
			else if (fabsf( z ) > 1.e0)
			{
				c = 1.e0/z;
				s = sqrtf( (1.e0 - c)*(1.e0 + c) );
			}
			else if (z == 1.e0)
			{
				c = 0.e0;
				s = 1.e0;
			}
			z = c*Delta[i] - s*Delta[j];
			Delta[j] = s*Delta[i] + c*Delta[j];
			Delta[i] = z;
		}
	}
 
 
	/*     Compute the amount of inconsistency in the constraint
	 *     equations.  Return a factor that allows for testing. */
	if (nrank < nc)
	{
		tn = snrm2( nc - nrank, &Delta[neq + nrank + 1], 1 );
	}
	else
	{
		tn = 0.e0;
	}
	/*     Store the value in delta(neq+1) for testing after returning.
	 *     Make that error test based on the size of delta(neq+1). */
	Delta[neq + 1] = tn;
	/*     Apply weights used for step.
	 *     They were also input in DELTA(:). */
	for (j = 1; j <= neq; j++)
	{
		Delta[j] *= sqrtf( Wt[j] );
	}
	return;
#undef	A
} /* end of function */