tms2.c

求矩阵奇异分解svd算法

               for(i=0; i<n; i++) z[i] = z4[i];
             }
         }
        }
  daxpy(n,eps,y,1,z,1);
  return;
} 
#include <stdio.h>
#include <math.h>
#include "tmsc.h"


void pmul(int, float *, float *, float *, float *,
          float *, float, int , float);
void dgemv2(int transa, int m, int n, float alpha,
        float **a, int ira, int ica, float *x,
        float beta, float *y);
void daxpy(int, float, float *, int, float *, int);
float ddot(int, float *, int, float *, int);
void dscal(int, float , float *, int);

void porth(int p, int f, int n, float **x, float *tmp,
           float *tmp1, float *tmp2, float *tmp3,
           float *tmp4, float e, int degree, float alpha)
/* 
     C translation of blas2  Gram-Schmidt orthogonalization procedure 
     for polynomial accelerated trace minimization (job=2)
     the n by p matrix z stored in columns f+1 to f+p of
     array x is reorthogonalized w.r.t. the first f columns of
     array x.  the resulting matrix z contains the desired P-orthogonal
     columns, where P is a polynomial in the matrix b from tsvd2.
     (based on orthog from Rutishauser)  */

{
int fpp,k,km1,j,i;

float s,t;

  if(!p) return;
  fpp = f+p;

  for(k=f; k<fpp; k++) {
     km1 = k-1;
L10:
    t = 0.0;
    if(km1<0) goto L25;
    pmul(n,x[k],tmp1,tmp2,tmp3,tmp4,e,degree,alpha);

  if(km1>0) {
    dgemv2(TRANSP,n,km1+1,1.0,x,0,0,tmp1,0.0,tmp);
    t += ddot(km1+1,tmp,1,tmp,1);
  }
  else {
     tmp[0] = ddot(n,x[0],1,tmp1,1);
     t += tmp[0]*tmp[0];
  }

  if(km1>0) dgemv2(NTRANSP,n,km1+1,-1.0,x,0,0,tmp,1.0,x[k]);
  else daxpy(n,-tmp[0],x[0],1,x[k],1);
  if(p==0) goto L50;
L25:
  pmul(n,x[k],tmp1,tmp2,tmp3,tmp4,e,degree,alpha);
  s = ddot(n,x[k],1,tmp1,1);
  t += s;
  if(s>t/100.0) goto L40;
  goto L10;
L40:
  s = sqrt(s);
  if(s!=0.0) s = 1.0/s;
  dscal(n,s,x[k],1);
L50:
  j=0;
  }
   return;
}
#include <stdio.h>

extern int ncol, nrow;
extern int *pointr, *rowind;
extern float *value ;

/************************************************************** 
 * multiplication of matrix B by a vector x, where            *
 *							      *
 * B =  A'A, where A is nrow by ncol (nrow >> ncol)           *
 * Hence, B is of order n:=ncol                               *
 * y stores product vector                                    *
 **************************************************************/ 
void opa(float *x, float *y)

{
   int i, j,end;

   for (i = 0; i < nrow; i++) y[i] = ZERO;

/* multiply by sparse C */
   for (i = 0; i < ncol; i++) {
      end = pointr[i+1];
      for (j = pointr[i]; j < end; j++)
	y[rowind[j]] += value[j] * x[i];
   }

   return;
}
#include <stdio.h>
#include <math.h>

void dtrsm(char side, char uplo, int transa, char diag, int m,
           int n, float alpha, float **a, float **b)

/*  dtrsm  solves one of the matrix equations
 *
 *     op( a )*x = alpha*b,   or   x*op( a ) = alpha*b,
 *
 *  where alpha is a scalar, x and b are m by n matrices, a is a unit, or
 *  non-unit,  upper or lower triangular matrix  and  op( a )  is one  of
 *
 *     op( a ) = a   or   op( a ) = a'.
 *
 *  the matrix x is overwritten on b.
 *
 *  parameters
 *  ==========
 *
 *  side   - character*1.
 *           on entry, side specifies whether op( a ) appears on the left
 *           or right of x as follows:
 *
 *              side = 'l' or 'l'   op( a )*x = alpha*b.
 *
 *              side = 'r' or 'r'   x*op( a ) = alpha*b.
 *
 *           unchanged on exit.
 *
 *  uplo   - character*1.
 *           on entry, uplo specifies whether the matrix a is an upper or
 *           lower triangular matrix as follows:
 *
 *              uplo = 'u' or 'u'   a is an upper triangular matrix.
 *
 *              uplo = 'l' or 'l'   a is a lower triangular matrix.
 *
 *           unchanged on exit.
 *
 *  transa - int.
 *           on entry, transa specifies the form of op( a ) to be used in
 *           the matrix multiplication as follows:
 *
 *              transa = 0    op( a ) = a.
 *
 *              transa = 1    op( a ) = a'.
 *
 *
 *           unchanged on exit.
 *
 *  diag   - character*1.
 *           on entry, diag specifies whether or not a is unit triangular
 *           as follows:
 *
 *              diag = 'u' or 'u'   a is assumed to be unit triangular.
 *
 *              diag = 'n' or 'n'   a is not assumed to be unit
 *                                  triangular.
 *           unchanged on exit.
 *
 *  m      - integer.
 *           on entry, m specifies the number of rows of b. m must be at
 *           least zero.
 *           unchanged on exit.
 *
 *  n      - integer.
 *           on entry, n specifies the number of columns of b.  n must be
 *           at least zero.
 *           unchanged on exit.
 *
 *  alpha  - float precision.
 *           on entry,  alpha specifies the scalar  alpha. when  alpha is
 *           zero then  a is not referenced and  b need not be set before
 *           entry.
 *           unchanged on exit.
 *
 *  a      - float precision array of dimension ( lda, k ), where k is m
 *           when  side = 'l' or 'l'  and is  n  when  side = 'r' or 'r'.
 *           before entry  with  uplo = 'u' or 'u',  the  leading  k by k
 *           upper triangular part of the array  a must contain the upper
 *           triangular matrix  and the strictly lower triangular part of
 *           a is not referenced.
 *           before entry  with  uplo = 'l' or 'l',  the  leading  k by k
 *           lower triangular part of the array  a must contain the lower
 *           triangular matrix  and the strictly upper triangular part of
 *           a is not referenced.
 *           note that when  diag = 'u' or 'u',  the diagonal elements of
 *           a  are not referenced either,  but are assumed to be  unity.
 *           unchanged on exit.
 *
 *  lda    - integer.
 *           on entry, lda specifies the first dimension of a as declared
 *           in the calling (sub) program.  when  side = 'l' or 'l'  then
 *           lda  must be at least  max( 1, m ),  when  side = 'r' or 'r'
 *           then lda must be at least max( 1, n ).
 *           unchanged on exit.
 *
 *  b      - float precision array of dimension ( ldb, n ).
 *           before entry,  the leading  m by n part of the array  b must
 *           contain  the  right-hand  side  matrix  b,  and  on exit  is
 *           overwritten by the solution matrix  x.
 *
 *  ldb    - integer.
 *           on entry, ldb specifies the first dimension of b as declared
 *           in  the  calling  (sub)  program.   ldb  must  be  at  least
 *           max( 1, m ).
 *           unchanged on exit.
 *
 *
 *  C translation of level 3 blas routine.
 *  -- written on 8-february-1989.
 *     jack dongarra, argonne national laboratory.
 *     iain duff, aere harwell.
 *     jeremy du croz, numerical algorithms group ltd.
 *     sven hammarling, numerical algorithms group ltd.   */
 
{
int i,j,k,lside,nounit,upper,nrowa;
float temp;

  if(side=='L') { 
    nrowa = m;
    lside=1;
  }
  else {
    nrowa = n;
    lside = 0;
  }
  if(diag=='N') nounit = 1; 
  else nounit = 0;
  if(uplo == 'U') upper = 1;
  else upper = 0;

  /*info = 0;
  if((!lside) && (!lsame(side,'R')) info =1; */
  
  
  if(alpha==ZERO) {
  for(j=0; j<n; j++)
     for(i=0; i<m; i++)
        b[j][i] = ZERO;
  return;
  }
if(lside) {
  if(!transa) {
    if(upper) {
      for(j=0; j<n; j++) {
         if(alpha!=ONE) for(i=0; i<m; i++) b[j][i] *= alpha;
         for(k=m-1; k>=0; k--) {
            if(b[j][k]!=ZERO) {
              if(nounit) b[j][k] /= a[k][k];
              for(i = 0; i<k; i++) b[j][i] -= b[j][k]*a[k][i];
            }
         }
      }
     }
    else {
     for(j=0; j<n; j++) {
        if(alpha!=ONE) for(i=0; i<m; i++) b[j][i] *= alpha;
        for(k=0; k<m; k++) {
             if(b[j][k]!=ZERO) {
               if(nounit) b[j][k] /= a[k][k];
               for(i=k+1; i<m; i++) b[j][i] -= b[j][k]*a[k][i];
             }
        }
     }
    }
 }
 else {
/* b = alpha*inv(a') *b */

    if(upper) {
      for(j=0; j<n; j++) {
         for(i=0; i<m; i++) {
            temp = alpha*b[j][i];
            for(k=0; k<i; k++) temp -= a[i][k]*b[j][k];
            if(nounit) temp /= a[i][i];
            b[j][i] = temp;
         }
      }
    }
    else {
       for(j=0; j<n; j++) {
          for(i=m-1; i>=0; i--) {
              temp = alpha*b[j][i];
              for(k=i+1; k<m; k++) temp -= a[i][k]*b[j][k];
              if(nounit) temp /= a[i][i];
              b[j][i] = temp;
          }
       }
    }
 }
}
else {
/* b = alpha*b*inv(a) */
     if(!transa) {
       if(upper) {
         for(j=0; j<n; j++) { 
            if(alpha!=ONE) for(i=0; i<m; i++) b[j][i] *= alpha;
            for(k=0; k<j; k++) {
               if(a[j][k]!=ZERO) {
                 for(i=0; i<m; i++) b[j][i] -= a[j][k]*b[k][i];
               }
            }
            if(nounit) {
              temp = ONE/a[j][j];
              for(i=0; i<m; i++) b[j][i] *= temp;
            }
         }
        }
      else {
        for(j=n-1; j>=0; j--) {
           if(alpha!=ONE) for(i=0; i<m; i++) b[j][i] *= alpha;
           for(k=j+1; k<n; k++) {
              if(a[j][k] !=ZERO) 
                for(i=0; i<m; i++) b[j][i] -= a[j][k]*b[k][i];
            }
           if(nounit) {
             temp = ONE/a[j][j];
             for(i=0; i<m; i++) b[j][i] *= temp;
           }
        }
      }
     }
   else {
/* form b = alpha*b*inv[a']. */
   if(upper) {
    for(k=n-1; k>=0; k--) {
       if(nounit) {
         temp = ONE/a[k][k];
         for(i=0; i<m; i++) b[k][i] *= temp;
       }
       for(j=0; j<k; j++) {
          if(a[k][j]!=ZERO) {
            temp = a[k][j];
            for(i=0; i<m; i++) b[j][i] -=temp*b[k][i];
          }
       }
       if(alpha !=ONE) for(i=0; i<m; i++) b[k][i] *= alpha;
    }
   } 
   else {
      for(k=0; k<n; k++) {
         if(nounit) {
           temp = ONE/a[k][k];
           for(i=0; i<m; i++) b[k][i] *= temp;
         }
         for(j=k+1; j<n; j++) {
            if(a[k][j]!=ZERO) {
              temp = a[k][j];
              for(i=0; i<m; i++) b[j][i] -= temp*b[k][i];
            }
         }
         if(alpha!=ONE) for(i=0; i<m; i++) b[k][i] *= alpha;
     }
   }
 }
}
}

#include <math.h>
#include "tmsc.h"
float fsign(float, float);
float pythag(float, float);

/***********************************************************************
 *                                                                     *
 *				tql2()   			       *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   tql2() is a translation of a Fortran version of the Algol
   procedure TQL2, Num. Math. 11, 293-306(1968) by Dowdler, Martin, 
   Reinsch and Wilkinson.
   Handbook for Auto. Comp., vol.II-Linear Algebra, 227-240(1971).  

   This function finds the eigenvalues and eigenvectors of a symmetric
   tridiagonal matrix by the QL method.


   Arguments
   ---------

   (input)                                                             
   offset the index of the leading element  of the input(full) matrix
          to be factored.
   n      order of the symmetric tridiagonal matrix           
   d      contains the diagonal elements of the input matrix        
   e      contains the subdiagonal elements of the input matrix in its
            first n-1 positions.
   z      contains the identity matrix				    
                                                                   
   (output)                                                       
   d      contains the eigenvalues in ascending order.  if an error
            exit is made, the eigenvalues are correct but unordered for
            for indices 0,1,...,ierr.				   
   e      has been destroyed.					  
   z      contains orthonormal eigenvectors of the symmetric   
            tridiagonal (or full) matrix.  if an error exit is made,
            z contains the eigenvectors associated with the stored 
          eigenvalues.					
   ierr   set to zero for normal return, j if the j-th eigenvalue has
            not been determined after 30 iterations.		    
	  (return value)


   Functions used
   --------------
   UTILITY	fsign
   MISC		pythag

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

int tql2(int offset, int n, float *d, float *e, float **z)

{
   int j, last, l, l1, l2, m, i, k, iteration;
   float tst1, tst2, g, r, s, s2, c, c2, c3, p, f, h, el1, dl1;
   if (n == 1) return(0);
   f = ZERO;
   last = n - 1;
   tst1 = ZERO;
   e[last] = ZERO;

   for (l = offset; l < n; l++) {
      iteration = 0;
      h = fabs(d[l]) + fabs(e[l]);
      if (tst1 < h) tst1 = h;

      /* look for small sub-diagonal element */
      for (m = l; m < n; m++) {
	 tst2 = tst1 + fabs(e[m]);
	 if (tst2 == tst1) break;
      }
      if (m != l) {
	 while (iteration < 30) {