tms.cc

矩阵奇异分解(svd)最新c++版本

  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;
}


void tms::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;
     }
   }
 }
}
}


void tms::pmul(long n, double *y, double *z, double *z2, double *z3,
          double *z4, double e, long degree, double alpha)

/* Compute the multiplication z=P(B)*y (and leave y untouched).
   P(B) is constructed from chebyshev polynomials. The input
   parameter degree specifies the polynomial degree to be used.
   Polynomials are constructed on interval (-e,e) as in Rutishauser's
   ritzit code                                                    */
{

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) myopb(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);
          myopb(n,z,z3,alpha2);
          for(kk=1; kk<degree; kk++) {
             for(i=0; i<n; i++) {
                z4[i] = z3[i];
                z[i] = -z[i];
             }
             myopb(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;
}

void tms::myopb(long n, double *x, double *y, double shift)
{
/* multiplication of B = [(alpha*alpha-shift)*I - A'A] by a vector x , where

        where A is nrow by ncol (nrow>>ncol), and alpha an upper bound
        for the largest singular value of A, and 
        shift is the approximate singular value of A.
        (y stores product vector)                                          */

  long i,j,last;

  for(i=0; i<n; i++)
     y[i]=(alpha*alpha-shift)*x[i];

  for(i=0; i<nrow; i++) zz[i] = 0.0;

  for(i=0; i<ncol; i++) {
     last = pointr[i+1];
     for(j=pointr[i]; j<last; j++)
        zz[rowind[j]] += value[j]*x[i];
  }

  for(i=0; i<ncol; i++) {
     for(j=pointr[i]; j<pointr[i+1]; j++)
        y[i] -= value[j]*zz[rowind[j]];
  }
  return;
}