bls2.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

   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
            leading m by k part of a must contain the matrix A. Otherwise,
	    the leading k by m part of a must contain  the matrix A.

   b        matrix B as a 2-dimensional array.  When transb = NTRANSP, the
            leading k by n part of a must contain the matrix B. Otherwise,
	    the leading n by k part of a 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 2-dimensional array.  On entry, the leading
	    m by n part of c must contain the matrix C, except when
	    beta = 0.  In that case, c need not be set on entry. 
	    On exit, c is overwritten by the m by n matrix
	    (alpha * op(A) * op(B) + beta * C).

 ***********************************************************************/
void dgemm2(int transa, int transb, int m, int n, int k, 
            float alpha, float **a, float **b, float beta, float **c)
{
   int info;
   int i, j, l, nrowa, ncola, nrowb, ncolb;
   float temp, *atemp;

   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 %1d %s\n",
      "*** ON ENTRY TO DGEMM2, 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;
   }
   if (!m || !n || ((alpha == ZERO || !k) && beta == ONE))
      return;

   if (alpha == ZERO) {
      if (beta == ZERO) 
         for (i = 0; i < m; i++)
            for (j = 0; j < n; j++) c[i][j] = ZERO;

      else if (beta != ONE)
         for (i = 0; i < m; i++)
            for (j = 0; j < n; j++) c[i][j] *= beta;
      return;
   }

   if (beta == ZERO)
      for (i = 0; i < m; i++)
         for (j = 0; j < n; j++) c[i][j] = ZERO;

   else if (beta != ONE)
      for (i = 0; i < m; i++)
         for (j = 0; j < n; j++) c[i][j] *= beta;

   if (!transb) { 

      switch(transa) {

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

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

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

	 /* form C := alpha * A' * B' + beta * C */
	 case TRANSP:   for(i = 0; i < ncola; i++) {
			   for (l = 0; l < nrowb; l++) {
      	       		      temp = ZERO;
	 		      for(j = 0; j < nrowa; j++) 
				 temp += a[j][i] * b[l][j];
                              c[i][l] += 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 float precision.


   Parameters
   ----------

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

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

float enorm(int n, float *x)

{
   float norm2, agiant, floatn, s1, s2, s3, xabs, x1max, x3max;
   int i;

   s1     = ZERO;
   s2     = ZERO;
   s3     = ZERO;
   x1max  = ZERO;
   x3max  = ZERO;
   floatn = (float)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-triangular.

   p       number of block partitions used in blklan1.c.

		    s = s[0]               r= r[0]
			----                  ----
			s[1]                  r[1]
			----                  ----
			  .                     .
                          .                     .
			  .                     .
			s[p-1]                r[p-2]

   n       dimension of bigs

   (output)
   bigs    The 2-dimensional array bigs will contain this matrix in the 
	   following band matrix format:

           EXAMPLE WITH 4 SUP. DIAGONALS:
  
             TRANSPOSE OF   [0 0 0 0--1ST SUP. DIAGONAL----]
             ------------   [0 0 0 ---2ND SUP. DIAGONAL----]
                            [0 0 -----3RD SUP. DIAGONAL----]
                            [0 -------4TH SUP. DIAGONAL----]
                            [-------- MAIN DIAGONAL--------]
  
    Note that the super-diagonals and main diagonal of S are stored
    in COLUMNS of bigs (bigs is n by n).  

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

void formbigs(int n, int ndp, float **r, float **s, float **bigs)

{
   int p, i, j, k, kk, m, row, col;

   p = n / ndp;

   /* load main diagonal of bigs */
   j = 0;
   for (i = 0; i < p; i++) 
      for (k = 0; k < ndp; k++) { 
	 bigs[j][ndp] = s[j][k];
	 j++;
      }

   /* load super-diagonals of bigs (from top to bottom) */
   for (i = 0; i < ndp; i++) {

      /* pad zeros at start of a column */
      for (kk = 0; kk < ndp - i; kk++) 
	 bigs[kk][i] = ZERO;

      /* load first row of bigs with main diagonals of r[j] */
      if (i == 0) {
	 j = 0;
         for (m = 0; m < p - 1; m++) 
            for (k = 0; k < ndp; k++) 
	       bigs[kk++][0] = r[j++][k];
      }
      else {
	 m = 0;
         for (j = 0; j < p; j++) {
	    row = m;
	    col = ndp - i;
	    /* load elements form s[j] submatrices */
	    while (col < ndp)
	       bigs[kk++][i] = s[row++][col++];

	    /* load elements form r[j] submatrices */
            if (j < p - 1) {
	       col = i;
	       row = m;
	       while (col < ndp)
	         bigs[kk++][i] = r[row++][col++];
            }
            m += ndp;
         }
      }
   }
}
#include <stdio.h>
#define		ZERO	 0.0
#define		ONE	 1.0

int imax(int, int);

/***********************************************************************
 *                                                                     *
 *                          dsbmv()                                    *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   The function performs the matrix-vector operation

	  y := alpha * A * y + beta * y,

   where alpha and beta are scalars, x and y are n-element vectors and
   A is an n by n symmetric band matrix, with k super-diagonals.

   Parameters
   ----------

   n         number of rows of matrix A; n must be at least 0.  Unchanged
	     upon exit.

   k         number of super-diagonals of matrix A

   a         2-dimensional array whose leading n by (k + 1) part must 
	     contain the upper triangular band part of the symmetric matrix,
	     supplied row by row, with the leading diagonal of the matrix 
	     in column (k + 1) of the array, the first super-diagonal 
	     starting at position 2 in column k, and so on.
	     The top left k by k triangle of the array A is not referenced.

   x         linear array of dimension of at least n.  Before entry,
	     x must contain the n elements of vector x.  Unchanged on exit.

   y         linear array of dimension of at least n.  Before entry,
	     y must contain the n elements of vector y.  On exit, y is
	     overwritten by the updated vector y.


   Functions called
   --------------

   MISC      imax  

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

void dsbmv(int n, int k, float alpha, float **a,