s_bls2.cc

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

            ns = nc / nb;
         }
      }

      /*------------------------------------------------------------------*
       *                  PHASE 5 (back transformation)                   *
       *------------------------------------------------------------------*/

      if (nbold > 1) {
         dgemm2(NTRANSP, NTRANSP, nn, nn, nn, ONE, qqp, ppp, ZERO, ztempp);
         dgemm2(NTRANSP, NTRANSP, nb + k0, n, nn, ONE, ztempp, vvp, ZERO, uvtmpp);
	 ptr1 = uvtmp;
	 length = (nb + k0) * n;
	 for (i = 0; i < length; i++) v[i] = *ptr1++;
      }
      else dgemv(TRANSP, nn, n, ONE, vvp, qq, ZERO, v);
      if (k0) {
         final = iconv + k0;
         ii = 0;
         for (jj = iconv; jj < final; jj++) {
	    ptr1 =  &v[n * index[ii++]];
	    ptr2 = &v0[n * jj];
            for (i = 0; i < n; i++) *ptr2++ = *ptr1++;
         }
         iconv = final;
         if (k <= 0) {
	    iter -= 1;
	    cont = FALSE;
         }
         if (cont) {
            kk = 0;
            final = nb + k0;
            for (jj = 0; jj < final; jj++) {
	       flag = FALSE;
	       ll = 0;
	       while (ll < k0 && !flag) {
	          if (index[ll] != jj) ll++;
	          else flag = TRUE;
	       }
	       if (!flag) {
	          ptr1 = &vv[n * kk];
	          ptr2 =  &v[n * jj];
                  for (i = 0; i < n; i++) *ptr1++ = *ptr2++;
	          kk++;
	       }
            }
	 }
      }
      else { 
         ptr1 =  v;
         for (jj = 0; jj < nb; jj++)
            for (i = 0; i < n; i++) vvp[jj][i] = *ptr1++; 
      }
      if (!cont) break;
   }
   ptr1 =  v;
   ptr2 = v0;
   for (j = 0; j < iconv; j++) 
      for (i = 0; i < n; i++) *ptr1++ = *ptr2++;
   *iko = ik - k;
   *ibo  = nb;
   *ico = nc;

   free(workptr1);
   free(workptr2);
   free(index);
   return(0);
}


/***********************************************************************
 *                                                                     *
 *                        polong2()                                     *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   This function is a single-vector Lanczos tridiagonalization procedure
   (degenerate case of block size = 1 for block2.c) which is used when
   normal deflation reduces the current block size to 1.

   The function returns DONE when the number of remaining triplets to be 
   approximated is less than or equal to zero.  Otherwise, it returns
   status CONTINUE.
   

   Arguments
   ---------

   (input)
   k        current # of desired triplets
   m        row dimension of the sparse matrix A whose SVD is sought
   n        column dimension of the sparse matrix A whose SVD is sought
   tol     user-specified tolerance for approximate singular triplets
   wp      work space

   (output)
   sing    linear array containing the iko approximate singular values
   res     linear array containing the iko residuals of the approximate 
	   singular triplets
   alpha   diagonal elements of symmetric tridiagonal matrix from the 
	   inner recursion
   beta    off-diagonal elements of symmetric tridiagonal matrix from the 
	   inner recursion

   External parameters
   -------------------

   Defined and documented in bls2.h


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

   BLAS         ddot, enorm, daxpy, orthg
   USER         opb

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

long s_bls2::polong2(long *k, long m, long n, 
            double *v, double *res, double *eig, double tol)

{
   double *ptr, dum, znm;
   long i, j, jj, convpj;

   for (j = 0; j < n; j++) p[j] = vv[j];
   for (j = 0; j < n; j++) z[j] = ZERO;

   for (j = 0; j < nn; j++) {
      opb(m, n, p, q);
      daxpy(n, ONE, q, 1, z, 1);
      alpha[j] = ddot(n, p, 1, z, 1);
      if (j == nn - 1) return(CONTINUE);

      /* compute Z[j] := Z[j] - alpha[j] * P */
      if (alpha[j] != ZERO) {
	 daxpy(n, -alpha[j], p, 1, z, 1);
	 if (!j) {
	    if ((znm = enorm(n, z)) <= tol) {
	       ptr = &v[iconv * n];
               for (i = 0; i < n; i++) *ptr++ = p[i];
	       eig[iconv] = alpha[0];
	       res[iconv] = znm;
               *k -= 1;
	       iter -= 1;
	       return(DONE);
            }
	 }

	 /* orthogonalize Z w.r.t. converged right S-vectors and previous 
	    VV's */
         convpj = iconv + j;
	 ptr = v0;
	 for (jj = 0; jj < iconv; jj++)
	    for (i = 0; i < n; i++)
	       uvtmpp[jj][i] = *ptr++;
	 ptr = vv;
	 for (jj = iconv; jj <= convpj; jj++)
	    for (i = 0; i < n; i++)
	       uvtmpp[jj][i] = *ptr++;
         if (convpj) {
	    ptr = uvtmpp[convpj + 1];
	    for (i = 0; i < n; i++) *ptr++ = z[i];
	    orthg(1, convpj + 1, n, yp, uvtmpp, ztemp);
	    ptr = uvtmpp[convpj + 1];
	    for (i = 0; i < n; i++) z[i] = *ptr++;
	 }
	 else {
	    dum = -ddot(n, uvtmp, 1, z, 1);
	    daxpy(n, dum, uvtmp, 1, z, 1);
	 }

	 /* compute beta[j] */
	 beta[j] = enorm(n,z);
	 if (beta[j] != ZERO) {

            /* compute P[j+1] := Z[j] / beta[j] */
	    ptr = vvp[j + 1];
	    for (i = 0; i < n; i++) {
	       t[i] = p[i];
	       p[i] = z[i] / beta[j];
	       *ptr++ = p[i];
	       z[i] = -beta[j] * t[i];
	    }
         }
	 else return(CONTINUE);
      }
      else return(CONTINUE);
   }
}


/***********************************************************************
 *                                                                     *
 *                        block2()                                     *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   This function implements the first two phases of the hybrid block
   Lanczos procedure.  In the first phase, which is also known as the
   block Lanczos outer iteration, a symmetric block tridiagonal matrix S
   is formed.  The eigenvalues of the matrix S approximate those 
   of matrix B, where B = A'A and A is the original sparse matrix.  
   Total (or complete) re-orthogonalization is used.

   In the second phase, single-vector Lanczos tridiagonalization is used
   to reduce (preserving eigenvalues of S) the block matrix S to a 
   symmetric tridiagonal matrix T. 


   Arguments
   ---------

   (input)
   w, wp     work space
   sp        diagonal blocks (symmetric submatrices) of the
		symmetric block tridiagonal matrix S 
   rp        super-diagonal blocks (upper-triangular submatrices)
		of the symmetric block tridiagonal matrix S 
   bigsp     symmetric block tridiagonal matrix S
   m         row dimension of sparse matrix A
   n         column dimension of sparse matrix A
   nb        current block size
   ns        number of blocks in current iteration
   irand     seed for random number generator
   

   (output - globally defined)
   alpha     diagonal elements of symmetric tridiagonal matrix T (reduced 
                from matrix S)
   beta      off-diagonal elements of symmetric tridiagonal matrix T 
		(reduced from matrix S)
   tres      residuals of approximate eigenvalues determined from
	     a previous set of block Lanczos outer iterations.
   ppp       array of eigenvectors of S


   External parameters
   -------------------

   Defined and documented in bls2.h


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

   BLAS         ddot, daxpy, enorm, dgemm2, orthg, dsbmv
   USER         opm
   MISC         random
   BLS2         formbigs

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

void s_bls2::block2(double **sp, double **rp, double **bigsp, long m, long n, 
            long nb, long ns, long *irand)

{
   long jinc, nk, nj, i, j, k, blkptr;
   double *ptr, dum, pnm;

   for (i = 0; i < nn; i++)
      for (j = 0; j < nb; j++) sp[i][j] = ZERO;
   /* ns (number of blocks) is assumed to be at least 2 */
   nk = nn - nb;
   for (i = 0; i < nk; i++)
      for (j = 0; j < nb; j++) rp[i][j] = ZERO;

   opm(m, n, nb, vvp, yp);

   dgemm2(NTRANSP, TRANSP, nb, nb, n, ONE, vvp, yp, ZERO, sp);

   blkptr = 0;
   for (j = 1; j < ns; j++) {

      dgemm2(TRANSP, NTRANSP, nb, n, nb, -ONE, &sp[blkptr],
            &vvp[blkptr], ONE, yp);
      if (j > 1)
         dgemm2(NTRANSP, NTRANSP, nb, n, nb, -ONE, &rp[blkptr - nb],
               &vvp[blkptr - nb], ONE, yp);

      if (j == 1 && iter > 1)
         for (i = 0; i < nb; i++) tres[i] = enorm(n, yp[i]);

      if (iconv) {
         nk = nb * j;
         nj = nk + iconv;

         ptr = vv;
         for (k = 0; k < nk; k++) 
            for (i = 0; i < n; i++) uvtmpp[k][i] = *ptr++; 
         ptr = v0;
         for (k = nk; k < nj; k++) 
            for (i = 0; i < n; i++) uvtmpp[k][i] = *ptr++; 
      }
      else {
         nj = nb * j;
         ptr = vv;
         for (k = 0; k < nj; k++) 
            for (i = 0; i < n; i++) uvtmpp[k][i] = *ptr++;
      }
      ptr = y;
      for (k = 0; k < nb; k++) 
         for (i = 0; i < n; i++) {
	    uvtmpp[nj + k][i] = *ptr;
	    *ptr++ = ZERO;
         }
      orthg(nb, nj, n, yp, uvtmpp, ztemp);

      for (k = 0; k < nb; k++) 
         for (i = k; i < nb; i++) rp[blkptr + k][i] = yp[i][k];

      jinc = blkptr + nb;
      ptr = vvp[jinc];
      for (k = nj; k < nj + nb; k++) 
         for (i = 0; i < n; i++) *ptr++ = uvtmpp[k][i]; 
  
      opm(m, n, nb, &vvp[jinc], yp);
      dgemm2(NTRANSP, TRANSP, nb, nb, n, ONE, &vvp[jinc], yp, ZERO, &sp[jinc]);
      blkptr += nb;
   }
   formbigs(nn, nb, rp, sp, bigsp);

   for (i = 0; i < nn; i++) p[i] = mrandom(irand);
   pnm = enorm(nn,p);

   ptr = pp;
   for (i = 0; i < nn; i++) {
      p[i] /= pnm;
      *ptr++ = p[i];
      z[i] = ZERO;
   }
   for (j = 0; j < nn; j++) {
      dsbmv(nn, nb, ONE, bigsp, p, ONE, z);
      alpha[j] = ddot(nn, p, 1, z, 1);
      if (j == nn - 1) break;

      /* compute Z[j] := Z[j] - alpha[j] * P */
      daxpy(nn, -alpha[j], p, 1, z, 1);

      /* orthogonalize Z w.r.t. previous PP's */
      for (k = 0; k <= j; k++)
         for (i = 0; i < nn; i++) uvtmpp[k][i] = ppp[k][i];
      if (j) {
         ptr = uvtmpp[j + 1];
         for (i = 0; i < nn; i++) *ptr++ = z[i];
	 orthg(1, j+1, nn, yp, uvtmpp, ztemp);
	 ptr = uvtmpp[j + 1];
         for (i = 0; i < nn; i++) z[i] = *ptr++;
      }
      else {
         dum = -ddot(nn, uvtmp, 1, z, 1);
	 daxpy(nn, dum, uvtmp, 1, z, 1);
      }
      beta[j] = enorm(nn,z);
      if (beta[j] != ZERO) {

	 /* compute P[j+1] := Z[j] / beta[j] */
         ptr = ppp[j + 1];
         for (i = 0; i < nn; i++) {
	    t[i] = p[i];
	    p[i] = z[i] /beta[j];
	    *ptr++ = p[i];
	    z[i] = -beta[j] * t[i];
	 }
      }
      else return;
   }
   return;
}