clinpack.c

来自「开放源码的编译器open watcom 1.6.0版的源代码」· C语言 代码 · 共 1,405 行 · 第 1/3 页

		 t = a[lda*k+l];
		 a[lda*k+l] = a[lda*k+k];
		 a[lda*k+k] = t; 
	    }

	    /* compute multipliers */

	    t = -ONE/a[lda*k+k];
	    dscal(n-(k+1),t,&a[lda*k+k+1],1);

	    /* row elimination with column indexing */

	    for (j = kp1; j < n; j++) {
		 t = a[lda*j+l];
		 if (l != k) {
		  a[lda*j+l] = a[lda*j+k];
		  a[lda*j+k] = t;
		 }
		 daxpy(n-(k+1),t,&a[lda*k+k+1],1,
		     &a[lda*j+k+1],1);
		} 
	   }
	 else { 
		     *info = k;
	 }
	} 
   }
   ipvt[n-1] = n-1;
   if (a[lda*(n-1)+(n-1)] == ZERO) *info = n-1;
}


/*----------------------*/ 

void dgesl(a,lda,n,ipvt,b,job)
int lda,n,ipvt[],job;
REAL a[],b[];

/* We would like to declare a[][lda], but c does not allow it.  In this
function, references to a[i][j] are written a[lda*i+j].  */

/*
     dgesl solves the double precision system
     a * x = b  or  trans(a) * x = b
     using the factors computed by dgeco or dgefa.

     on entry

	a       double precision[n][lda]
		the output from dgeco or dgefa.

	lda     integer
		the leading dimension of the array  a .

	n       integer
		the order of the matrix  a .

	ipvt    integer[n]
		the pivot vector from dgeco or dgefa.

	b       double precision[n]
		the right hand side vector.

	job     integer
		= 0         to solve  a*x = b ,
		= nonzero   to solve  trans(a)*x = b  where
			    trans(a)  is the transpose.

    on return

	b       the solution vector  x .

     error condition

	a division by zero will occur if the input factor contains a
	zero on the diagonal.  technically this indicates singularity
	but it is often caused by improper arguments or improper
	setting of lda .  it will not occur if the subroutines are
	called correctly and if dgeco has set rcond .gt. 0.0
	or dgefa has set info .eq. 0 .

     to compute  inverse(a) * c  where  c  is a matrix
     with  p  columns
	   dgeco(a,lda,n,ipvt,rcond,z)
	   if (!rcond is too small){
		for (j=0,j<p,j++)
		    dgesl(a,lda,n,ipvt,c[j][0],0);
	}

     linpack. this version dated 08/14/78 .
     cleve moler, university of new mexico, argonne national lab.

     functions

     blas daxpy,ddot
*/


{
/*     internal variables   */

   REAL ddot(),t;
   int k,kb,l,nm1;

   nm1 = n - 1;
   if (job == 0) {

	/* job = 0 , solve  a * x = b
	 first solve  l*y = b       */

	if (nm1 >= 1) {
	 for (k = 0; k < nm1; k++) {
	    l = ipvt[k];
	    t = b[l];
	    if (l != k){ 
		 b[l] = b[k];
		 b[k] = t;
	    }   
	    daxpy(n-(k+1),t,&a[lda*k+k+1],1,&b[k+1],1);
	 }
	} 

	/* now solve  u*x = y */

	for (kb = 0; kb < n; kb++) {
	  k = n - (kb + 1);
	  b[k] = b[k]/a[lda*k+k];
	  t = -b[k];
	  daxpy(k,t,&a[lda*k+0],1,&b[0],1);
	}
   }
   else { 

	/* job = nonzero, solve  trans(a) * x = b
	 first solve  trans(u)*y = b          */

	for (k = 0; k < n; k++) {
	 t = ddot(k,&a[lda*k+0],1,&b[0],1);
	 b[k] = (b[k] - t)/a[lda*k+k];
	}

	/* now solve trans(l)*x = y   */

	if (nm1 >= 1) {
	 for (kb = 1; kb < nm1; kb++) {
	    k = n - (kb+1);
	    b[k] = b[k] + ddot(n-(k+1),&a[lda*k+k+1],1,&b[k+1],1);
	    l = ipvt[k];
	    if (l != k) {
		 t = b[l];
		 b[l] = b[k];
		 b[k] = t;
	    }
	 }
	}
   }
}

/*----------------------*/ 

void daxpy(n,da,dx,incx,dy,incy)
/*
     constant times a vector plus a vector.
     jack dongarra, linpack, 3/11/78.
*/
REAL dx[],dy[],da;
int incx,incy,n;
{
   int i;
   int ix;
   int iy;
#ifdef UNROLL
   int m;
#endif
#if 0
   int mp1;
#endif

   if(n <= 0) return;
   if (da == ZERO) return;

   if(incx != 1 || incy != 1) {

	/* code for unequal increments or equal increments
	 not equal to 1                */

	ix = 1;
	iy = 1;
	if(incx < 0) ix = (-n+1)*incx + 1;
	if(incy < 0) iy = (-n+1)*incy + 1;
	for (i = 0;i < n; i++) {
	 dy[iy] = dy[iy] + da*dx[ix];
	 ix = ix + incx;
	 iy = iy + incy;
	}
	    return;
   }

   /* code for both increments equal to 1 */

#ifdef ROLL
   for (i = 0;i < n; i++) {
	dy[i] = dy[i] + da*dx[i];
   }
#endif
#ifdef UNROLL

   m = n % 4;
   if ( m != 0) {
	for (i = 0; i < m; i++) 
	 dy[i] = dy[i] + da*dx[i];
	if (n < 4) return;
   }
   for (i = m; i < n; i = i + 4) {
	dy[i]   = dy[i]   + da*dx[i];
	dy[i+1] = dy[i+1] + da*dx[i+1];
	dy[i+2] = dy[i+2] + da*dx[i+2];
	dy[i+3] = dy[i+3] + da*dx[i+3];
   }
#endif
}
   
/*----------------------*/ 

REAL ddot(n,dx,incx,dy,incy)
/*
     forms the dot product of two vectors.
     jack dongarra, linpack, 3/11/78.
*/
REAL dx[],dy[];

int incx,incy,n;
{
   REAL dtemp;
   int i;
   int ix;
   int iy;
#ifdef UNROLL
   int m;
#endif
#if 0
   int mp1;
#endif

   dtemp = ZERO;

   if(n <= 0) return(ZERO);

   if(incx != 1 || incy != 1) {

	/* code for unequal increments or equal increments
	 not equal to 1               */

	ix = 0;
	iy = 0;
	if (incx < 0) ix = (-n+1)*incx;
	if (incy < 0) iy = (-n+1)*incy;
	for (i = 0;i < n; i++) {
	 dtemp = dtemp + dx[ix]*dy[iy];
	 ix = ix + incx;
	 iy = iy + incy;
	}
	return(dtemp);
   }

   /* code for both increments equal to 1 */

#ifdef ROLL
   for (i=0;i < n; i++)
	dtemp = dtemp + dx[i]*dy[i];
   return(dtemp);
#endif
#ifdef UNROLL

   m = n % 5;
   if (m != 0) {
	for (i = 0; i < m; i++)
	 dtemp = dtemp + dx[i]*dy[i];
	if (n < 5) return(dtemp);
   }
   for (i = m; i < n; i = i + 5) {
	dtemp = dtemp + dx[i]*dy[i] +
	dx[i+1]*dy[i+1] + dx[i+2]*dy[i+2] +
	dx[i+3]*dy[i+3] + dx[i+4]*dy[i+4];
   }
   return(dtemp);
#endif
}

/*----------------------*/ 
void dscal(n,da,dx,incx)

/*     scales a vector by a constant.
	jack dongarra, linpack, 3/11/78.
*/
REAL da,dx[];
int n, incx;
{
   int i;
#if 0
   int m;
#endif
#ifdef UNROLL
   int mp1;
#endif
   int nincx;

   if(n <= 0)return;
   if(incx != 1) {

	/* code for increment not equal to 1 */

	nincx = n*incx;
	for (i = 0; i < nincx; i = i + incx)
	 dx[i] = da*dx[i];
	return;
   }

   /* code for increment equal to 1 */

#ifdef ROLL
   for (i = 0; i < n; i++)
	dx[i] = da*dx[i];
#endif
#ifdef UNROLL

   m = n % 5;
   if (m != 0) {
	for (i = 0; i < m; i++)
	 dx[i] = da*dx[i];
	if (n < 5) return;
   }
   for (i = m; i < n; i = i + 5){
	dx[i] = da*dx[i];
	dx[i+1] = da*dx[i+1];
	dx[i+2] = da*dx[i+2];
	dx[i+3] = da*dx[i+3];
	dx[i+4] = da*dx[i+4];
   }
#endif

}

/*----------------------*/ 
int idamax(n,dx,incx)

/*
     finds the index of element having max. absolute value.
     jack dongarra, linpack, 3/11/78.
*/

REAL dx[];
int incx,n;
{
   REAL dmax;
   int i, ix, itemp;

   if( n < 1 ) return(-1);
   if(n ==1 ) return(0);
   if(incx != 1) {

	/* code for increment not equal to 1 */

	ix = 1;
	dmax = fabs((double)dx[0]);
	ix = ix + incx;
	for (i = 1; i < n; i++) {
	 if(fabs((double)dx[ix]) > dmax)  {
	    itemp = i;
	    dmax = fabs((double)dx[ix]);
	 }
	 ix = ix + incx;
	}
   }
   else {

	/* code for increment equal to 1 */

	itemp = 0;
	dmax = fabs((double)dx[0]);
	for (i = 1; i < n; i++) {
	 if(fabs((double)dx[i]) > dmax) {
	    itemp = i;
	    dmax = fabs((double)dx[i]);
	 }
	}
   }
   return (itemp);
}

/*----------------------*/ 
REAL epslon (x)
REAL x;
/*
     estimate unit roundoff in quantities of size x.
*/

{
   REAL a,b,c,eps;
/*
     this program should function properly on all systems
     satisfying the following two assumptions,
	1.  the base used in representing dfloating point
	    numbers is not a power of three.
	2.  the quantity  a  in statement 10 is represented to 
	    the accuracy used in dfloating point variables
	    that are stored in memory.
     the statement number 10 and the go to 10 are intended to
     force optimizing compilers to generate code satisfying 
     assumption 2.
     under these assumptions, it should be true that,
	    a  is not exactly equal to four-thirds,
	    b  has a zero for its last bit or digit,
	    c  is not exactly equal to one,
	    eps  measures the separation of 1.0 from
		 the next larger dfloating point number.
     the developers of eispack would appreciate being informed
     about any systems where these assumptions do not hold.

     *****************************************************************
     this routine is one of the auxiliary routines used by eispack iii
     to avoid machine dependencies.
     *****************************************************************

     this version dated 4/6/83.
*/

   a = 4.0e0/3.0e0;
   eps = ZERO;
   while (eps == ZERO) {
	b = a - ONE;
	c = b + b + b;
	eps = fabs((double)(c-ONE));
   }
   return(eps*fabs((double)x));
}
 
/*----------------------*/ 
void dmxpy (n1, y, n2, ldm, x, m)
REAL y[], x[], m[];
int n1, n2, ldm;

/* We would like to declare m[][ldm], but c does not allow it.  In this
function, references to m[i][j] are written m[ldm*i+j].  */

/*
   purpose:
     multiply matrix m times vector x and add the result to vector y.

   parameters:

     n1 integer, number of elements in vector y, and number of rows in
	 matrix m

     y double [n1], vector of length n1 to which is added 
	 the product m*x

     n2 integer, number of elements in vector x, and number of columns
	 in matrix m

     ldm integer, leading dimension of array m

     x double [n2], vector of length n2

     m double [ldm][n2], matrix of n1 rows and n2 columns

 ----------------------------------------------------------------------
*/
{