bls1.c

求矩阵奇异分解svd算法

   ----------

   (input)
   transa   TRANSP indicates op(A) = A' is to be used in the multiplication
	    NTRANSP indicates op(A) = A is to be used in the multiplication

   m        on entry, m specifies the number of rows of the matrix A.
	    m must be at least zero.  Unchanged upon exit.

   n        on entry, n specifies the number of columns of the matrix A.
	    n must be at least zero.  Unchanged upon exit.

   alpha    a scalar multiplier.  Unchanged upon exit.

   a        matrix A as a 2-dimensional array.  Before entry, the leading
	    m by n part of the array a must contain the matrix A.

   x        linear array of dimension of at least n if transa = NTRANSP
	    and at least m otherwise.

   beta     a scalar multiplier.  When beta is supplied as zero then y
	    need not be set on input.  Unchanged upon exit.

   y        linear array of dimension of at least m if transa = NTRANSP
	    and at leat n otherwise.  Before entry with beta nonzero,
	    the array y must contain the vector y.  On exit, y is 
	    overwritten by the updated vector y.


 ***********************************************************************/

void dgemv(long transa, long m, long n, 
           double alpha, double **a, double *x, double beta, double *y)

{
   long info, leny, i, j;
   double temp, *ptrtemp;

   info = 0;
   if      ( transa != TRANSP && transa != NTRANSP ) info = 1;
   else if ( m < 0 ) 				     info = 2;
   else if ( n < 0 )				     info = 3;

   if (info) {
      fprintf(stderr, "%s %1ld %s\n",
      "*** ON ENTRY TO DGEMV, PARAMETER NUMBER",info,"HAD AN ILLEGAL VALUE");
      exit(info);
   }

   if (transa) leny = n;
   else        leny = m;

   if (!m || !n || (alpha == ZERO && beta == ONE))
      return;

   ptrtemp = y; 

   /* form Y := beta * Y */
   if (beta == ZERO) 
      for (i = 0; i < leny; i++) *ptrtemp++ = ZERO;
   else if (beta != ONE) 
      for (i = 0; i < leny; i++) *ptrtemp++ *= beta;

   if (alpha == ZERO) return;

   switch(transa) {

      /* form Y := alpha * A * X + Y */
      case NTRANSP:  for(i = 0; i < m; i++) {
                        ptrtemp = *a++;
		        temp = ZERO;
		        for(j = 0; j < n; j++) 
			   temp += *ptrtemp++ * x[j];
			y[i] += alpha * temp;
		     }
		     break;
		     
      /* form Y := alpha * A' * X + Y */
      case TRANSP:   for(i = 0; i < m; i++) { 
                        ptrtemp = *a++;
			if (x[i] != ZERO) {
			   temp = alpha * x[i];
			   for(j = 0; j < n; j++)
			      y[j] += temp * (*ptrtemp++);
			}
		     }
		     break;
   }
}
#include <stdio.h>

#define		TRANSP   1
#define		NTRANSP  0
#define		ZERO	 0.0
#define		ONE	 1.0

/***********************************************************************
 *                                                                     *
 *                         dgemm()                                     *
 *                                                                     *
 * A C-translation of the level 3 BLAS routine DGEMM by Dongarra,      *
 * Duff, du Croz, and Hammarling (see LAPACK Users' Guide).            *
 * In this version, two of the three arrays which store the matrices   *
 * used in this matrix-matrix multiplication are accessed as linear    *
 * arrays.                                                             *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

   Description
   -----------

   dgemm() performs one of the matrix-matrix operations

	      C := alpha * op(A) * op(B) + beta * C,

   where op(X) = X or op(X) = X', alpha and beta are scalars, and A, B
   and C are matrices, with op(A) an m by k matrix, op(B) a k by n
   matrix and C an m by n matrix.

   Note that the arrays storing matrices B and C are linear arrays while
   the array of A is two-dimensional.


   Parameters
   ----------

   (input)
   transa   TRANSP indicates op(A) = A' is to be used in the multiplication
	    NTRANSP indicates op(A) = A is to be used in the multiplication

   transb   TRANSP indicates op(B) = B' is to be used in the multiplication
	    NTRANSP indicates op(B) = B is to be used in the multiplication

   m        on entry, m specifies the number of rows of the matrix op(A)
	    and of the matrix C.  m must be at least zero.  Unchanged
	    upon exit.

   n        on entry, n specifies the number of columns of the matrix op(B)
	    and of the matrix C.  n must be at least zero.  Unchanged
	    upon exit.

   k        on entry, k specifies the number of columns of the matrix op(A)
	    and the number of rows of the matrix B.  k must be at least 
	    zero.  Unchanged upon exit.

   alpha    a scalar multiplier

   a        matrix A as a 2-dimensional array.  When transa = NTRANSP, the
	    first k columns of the first m rows must contain the matrix A.
	    Otherwise, the first m columns of the first k rows must contain
	    the matrix A.

   b        matrix B as a linear array.  The leading (k * n) elements of
	    b must contain the matrix B.

   beta     a scalar multiplier.  When beta is supplied as zero then C
	    need not be set on input.

   c        matrix C as a linear array.  Before entry, the leading (m * n)
	    elements of c must contain the matrix C except when beta = 0.
	    In this case, c need not be set on entry.
	    On exit, c is overwritten by the (m * n) elements of matrix
	    (alpha * op(A) * op(B) + beta * C).

 ***********************************************************************/

void dgemm(long transa, long transb, long m, long n, long k, 
           double alpha, double **a, double *b, double beta, double *c)

{
   long info;
   long i, j, l, nrowa, ncola, nrowb, ncolb, nc;
   double temp, *atemp, *btemp1, *ptrtemp, *ctemp;

   info = 0;
   if      ( transa != TRANSP && transa != NTRANSP ) info = 1;
   else if ( transb != TRANSP && transb != NTRANSP ) info = 2;
   else if ( m < 0 ) 				     info = 3;
   else if ( n < 0 )				     info = 4;
   else if ( k < 0 )        			     info = 5;

   if (info) {
      fprintf(stderr, "%s %1ld %s\n",
      "*** ON ENTRY TO DGEMM, PARAMETER NUMBER",info,"HAD AN ILLEGAL VALUE");
      exit(info);
   }

   if (transa) {
      nrowa = k;
      ncola = m;
   }
   else { 
      nrowa = m;
      ncola = k;
   }
   if (transb) {
      nrowb = n;
      ncolb = k;
   }
   else {
      nrowb = k;
      ncolb = n;
   }
   nc = m * n;

   if (!m || !n || ((alpha == ZERO || !k) && beta == ONE))
      return;

   ctemp = c; 
   if (alpha == ZERO) {
      if (beta == ZERO)
         for (i = 0; i < nc; i++) *ctemp++ = ZERO;
      else if (beta != ONE)
         for (i = 0; i < nc; i++) *ctemp++ *= beta;
      return;
   }

   if (beta == ZERO) 
      for (i = 0; i < nc; i++) *ctemp++ = ZERO;
   else if (beta != ONE) 
      for (i = 0; i < nc; i++) *ctemp++ *= beta;

   if (!transb) { 

      switch(transa) {

	 /* form C := alpha * A * B + beta * C */
	 case NTRANSP:  ptrtemp = c;
		        for(l = 0; l < nrowa; l++) {
                           atemp = *a++;
      	       		   btemp1 = b;
      		  	   for(j = 0; j < ncola; j++) {
	 	     	      temp = *atemp * alpha;
	             	      ctemp = ptrtemp;
         	     	      for(i = 0; i < ncolb; i++) 
	    			 (*ctemp++) += temp * (*btemp1++);
	 	     	      atemp++;
      		  	   }
   	       		   ptrtemp = ctemp;
   	    	        }
			break;

	 /* form C := alpha * A' * B + beta * C */
	 case TRANSP:   ptrtemp = b;
	                for(l = 0; l < nrowa; l++) {
                           atemp = *a++;
      	       		   ctemp = c;
      		  	   for(j = 0; j < ncola; j++) {
	 	     	      temp = *atemp * alpha;
	             	      btemp1 = ptrtemp;
         	     	      for(i = 0; i < ncolb; i++) 
	    			 (*ctemp++) += temp * (*btemp1++);
	 	     	      atemp++;
      		  	   }
   	       		   ptrtemp = btemp1;
   	    		}
			break;
      }
   }
   else { 
      ctemp = c;

      switch(transa) {

	 /* form C := alpha * A * B' + beta * C */
	 case NTRANSP: for(l = 0; l < nrowa; l++) {
      	       		   btemp1 = b;
      		  	   for(j = 0; j < nrowb; j++) {
	 	     	      atemp = *a;
         	     	      for(i = 0; i < ncolb; i++) 
	    			 *ctemp += (*atemp++) * alpha * (*btemp1++);
	 	     	      ctemp++;
      		  	   }
			   a++;
   	    		}
			break;

	 /* form C := alpha * A' * B' + beta * C */
	 case TRANSP:   for(i = 0; i < ncola; i++) {
			   btemp1 = b;
			   for (l = 0; l < nrowb; l++) {
      	       		      temp = ZERO;
	 		      for(j = 0; j < nrowa; j++) 
			         temp += a[j][i] * (*btemp1++);
	    		      *ctemp++ += alpha * temp;
			   }
   	    		}
			break;
      }
   }
}
#include <math.h>

#define       ZERO       0.0
#define       ONE        1.0
#define       RDWARF     3.834e-20
#define       RGIANT     1.304e19

/***********************************************************************
 *                                                                     *
 *                        enorm()                                      *
 *  a C translation of the Fortran-77 version by Burton, Garbow,       *
 *  Hillstrom and More of Argonne National Laboratory.                 *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

   Description
   -----------

   given an n-vector x, this function calculates the Euclidean norm of x.
   The Euclidean norm is computed by accumulating the sum of squares in 
   three different sums.  The sums of squares for the small and large 
   components are scaled so that no overflows occur.  Non-destructive 
   underflows are permitted.  Underflows and overflows do not occur in the 
   computation of the unscaled sum of squares for the intermediate components. 
   The definitions of small, intermediate and large components depend on two
   constants, rdwarf and rgiant.  The restrictions on these constants are 
   that rdwarf**2 not underflow and rgiant**2 not overflow.  The constants
   given here are suitable for every known computer.
   The function returns the Euclidean norm of vector x in double precision.


   Parameters
   ----------

   n         number of elements in vector x
   x         linear array of vector x whose Euclidean norm is to be 
		calculated

 ***********************************************************************/

double enorm(long n, double *x)

{
   double norm2, agiant, floatn, s1, s2, s3, xabs, x1max, x3max;
   long i;

   s1     = ZERO;
   s2     = ZERO;
   s3     = ZERO;
   x1max  = ZERO;
   x3max  = ZERO;
   floatn = (double)n;
   agiant = RGIANT / floatn;

   for (i = 0; i < n; i++) {
      xabs = fabs(x[i]);
      /* summing components of vector that need no scaling */
      if (xabs > RDWARF && xabs < agiant)
	 s2 += xabs * xabs;
      else {
	 /* underflow... */
	 if (xabs <= RDWARF) {
            if (xabs > x3max) {
	       s3 = ONE + s3 * (x3max/xabs) * (x3max/xabs);
	       x3max = xabs;
	    }
	    else if (xabs != 0)
	       s3 += (xabs/x3max) * (xabs/x3max);
	 }
	 /* overflow... */
	 else {
	    /* summing large components of vector */
	    if (xabs <= x1max)
	       s1 += (xabs/x1max) * (xabs/x1max);
            else {
	       s1 = ONE + s1 * (x1max/xabs) * (x1max/xabs);
	       x1max = xabs;
	    }
	 }
      }
   }
   if (s1 != ZERO)
      norm2 = x1max * sqrt(s1 + (s2/x1max) / x1max);
   else if (s2 != ZERO) {
      if (s2 >= x3max)
	 norm2 = sqrt(s2 * (ONE + (x3max/s2) * (x3max*s3)));
      else 
	 norm2 = sqrt(x3max * ((s2/x3max) + (x3max*s3)));
   }
   else
      norm2 = x3max * sqrt(s3);
   return(norm2);
}
#define         ZERO     0.0

/***********************************************************************
 *                                                                     *
 *                        formbigs()                                   *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

   Description
   -----------

   This function forms the block upper-bidiagonal or the symmetric block
   tridiagonal matrix S from the block Lanczos algorithm in Phase 1 of 
   blklan1.c or blklan2.c, respectively. 


   Arguments
   ---------

   (input)
   r, s    submatrices from which the bidiagonal block matrix S
	   (Phase 1 of blklan1.c) is formed.
           The following data structure is assumed for the submatrices
	   s[j] and r[j], where j = 0, 1, ..., p-1.  For blklan1.c,
	   s[j] and r[j] are both upper-triangular.  For blklan2.c,
	   s[j] is dense and symmetric and r[j] is upper-tria