tms2.c

求矩阵奇异分解svd算法

double tmp2,alpha2,eps,tmp;
 
long i,kk;

  alpha2 = alpha*alpha;
  tmp2 = 1.0;
  eps = tmp2 + 1.0e-6;
  if(!degree) for(i=0; i<n; i++) z[i] = y[i];
  else if(degree==1) opb(n,y,z,alpha2);
     else if(degree>1) {
          tmp2 = 2.0/e;
          tmp = 5.1e-1*tmp2;
          for(i=0; i<n; i++)
             z[i] = y[i];
          dscal(n,tmp,z,1);
          opb(n,z,z3,alpha2);
          for(kk=1; kk<degree; kk++) {
             for(i=0; i<n; i++) {
                z4[i] = z3[i];
                z[i] = -z[i];
             }
             opb(n,z3,z2,alpha2);
             daxpy(n,tmp2,z2,1,z,1);
             if(kk != degree-1) {
               for(i=0; i<n; i++) z3[i] = z[i];
               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(long, double *, double *, double *, double *,
          double *, double, long , double);
void dgemv2(long transa, long m, long n, double alpha,
        double **a, long ira, long ica, double *x,
        double beta, double *y);
void daxpy(long, double, double *, long, double *, long);
double ddot(long, double *, long, double *, long);
void dscal(long, double , double *, long);

void porth(long p, long f, long n, double **x, double *tmp,
           double *tmp1, double *tmp2, double *tmp3,
           double *tmp4, double e, long degree, double 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)  */

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

double 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 long ncol, nrow;
extern long *pointr, *rowind;
extern double *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(double *x, double *y)

{
   long 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, long transa, char diag, long m,
           long n, double alpha, double **a, double **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 - long.
 *           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  - double 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      - double 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      - double 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.   */
 
{
long i,j,k,lside,nounit,upper,nrowa;
double 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"
double fsign(double, double);
double pythag(double, double);

/***********************************************************************
 *                                                                     *
 *				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)