bls.cc

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

// File: bls.cc -- implements class bls
// Author: Suvrit Sra 
// The code is taken from svdpack....my contrib. is just the c++ wrapper.
// Date: 14 nov, 2003

/*************************************************************************
 THE ORIGINAL SVDPAKC COPYRIGHT
                           (c) Copyright 1993
                        University of Tennessee
                          All Rights Reserved                          
 *************************************************************************/

#include "bls.h"

using namespace ssvd;
/***********************************************************************
 *								       *
 *		              validate()			       *
 *								       *
 ***********************************************************************/
/***********************************************************************

   Description
   -----------
   Function validates input parameters and returns error code (long)  

   Arguments 
   ---------
  (input)
  fp_out1   pointer to output file (part of object)
  nnzero    number of nonzero elements in input sparse matrix A
  nrow      row dimension of A
  ncol      column dimension of A
  maxsubsp  maximum dimension of Krylov subspace allowed
  blksize   initial block size
  nums      number of singular values desired
  tol       user-specified tolerance for approximate singular triplets

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

long bls::validate(FILE* fp, long nnzero, long nrow, long ncol,
	      long maxsubsp, long blksize, long nums, double tol)
{
  long error_index;
  error_index = 0;
  if (ncol > NCMAX || nnzero > NZMAX)           error_index = 1;
  else if (ncol > nrow)                         error_index = 2;
  else if (tol < ZERO)                          error_index = 3;
  else if (maxsubsp > ncol || maxsubsp <= 0)    error_index = 4;
  else if (blksize <= 1 || blksize > maxsubsp)  error_index = 5;
  else if (nums > maxsubsp || nums <= 0)        error_index = 6;
  else if (blksize > nums)                      error_index = 7;
  
  if (error_index)
    fprintf(fp, "%s\n", error[error_index]);
  return(error_index);
}


/***********************************************************************
 *                                                                     *
 *                        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 bls::formbigs(long n, long ndp, double **r, double **s, double **bigs)

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