Subversion Repositories slepc-dev

/*                      

   SLEPc singular value solver: "lanczos"

   Method: Golub-Kahan-Lanczos bidiagonalization

   Last update: Nov 2006

*/

#include "src/svd/svdimpl.h"                /*I "slepcsvd.h" I*/

#include "slepcblaslapack.h"

#undef __FUNCT__  

#define __FUNCT__ "SVDSetUp_LANCZOS"

PetscErrorCode SVDSetUp_LANCZOS(SVD svd)

{

  PetscErrorCode  ierr;

  PetscInt        M,N;

  PetscFunctionBegin;

  ierr = MatGetSize(svd->A,&M,&N);CHKERRQ(ierr);

  if (svd->ncv == PETSC_DECIDE)

    svd->ncv = PetscMin(PetscMin(M,N),PetscMax(2*svd->nsv,10));

  if (svd->max_it == PETSC_DECIDE)

    svd->max_it = PetscMax(PetscMin(M,N)/svd->ncv,100);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "cgs"

PetscErrorCode cgs(Vec v, int n, Vec *V)

{

  PetscErrorCode ierr;

  Vec            w;

  PetscScalar    *h;

  PetscFunctionBegin;

  if (n>0) {

    ierr = VecDuplicate(v,&w);CHKERRQ(ierr);

    ierr = PetscMalloc(sizeof(PetscScalar)*n,&h);CHKERRQ(ierr);

    ierr = VecMDot(v,n,V,h);CHKERRQ(ierr);

    ierr = VecSet(w,0.0);CHKERRQ(ierr);

    ierr = VecMAXPY(w,n,h,V);CHKERRQ(ierr);

    ierr = VecAXPY(v,-1.0,w);CHKERRQ(ierr);

    ierr = VecDestroy(w);CHKERRQ(ierr);

    ierr = PetscFree(h);CHKERRQ(ierr);

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "computeres"

PetscErrorCode computeres(SVD svd,PetscScalar sigma,Vec u,Vec v,PetscReal *norm1,PetscReal *norm2)

{

  PetscErrorCode ierr;

  Vec            x,y;

  PetscFunctionBegin;

  ierr = VecDuplicate(u,&x);CHKERRQ(ierr);

  ierr = MatMult(svd->A,v,x);CHKERRQ(ierr);

  ierr = VecAXPY(x,-sigma,u);CHKERRQ(ierr);

  ierr = VecNorm(x,NORM_2,norm1);CHKERRQ(ierr);

  ierr = VecDuplicate(v,&y);CHKERRQ(ierr);

  if (svd->AT) {

    ierr = MatMult(svd->AT,u,y);CHKERRQ(ierr);

  } else {

    ierr = MatMultTranspose(svd->A,u,y);CHKERRQ(ierr);

  }

  ierr = VecAXPY(y,-sigma,v);CHKERRQ(ierr);

  ierr = VecNorm(y,NORM_2,norm2);CHKERRQ(ierr);

  ierr = VecDestroy(x);CHKERRQ(ierr);

  ierr = VecDestroy(y);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "SVDSolve_LANCZOS"

PetscErrorCode SVDSolve_LANCZOS(SVD svd)

{

  PetscErrorCode ierr;

  PetscReal      *alpha,*beta,norm1,norm2,*work;

  PetscScalar    *Q,*PT;

  int            i,j,k,l,n,zero=0,info;

  Vec            *V,*U;

  PetscFunctionBegin;

  ierr = PetscMalloc(sizeof(PetscReal)*svd->n,&alpha);CHKERRQ(ierr);

  ierr = PetscMalloc(sizeof(PetscReal)*svd->n,&beta);CHKERRQ(ierr);

  ierr = PetscMalloc(sizeof(PetscScalar)*svd->n*svd->n,&Q);CHKERRQ(ierr);

  ierr = PetscMalloc(sizeof(PetscScalar)*svd->n*svd->n,&PT);CHKERRQ(ierr);

  ierr = PetscMalloc(sizeof(PetscReal)*4*svd->n,&work);CHKERRQ(ierr);

  ierr = VecDuplicateVecs(svd->V[0],svd->n+1,&V);CHKERRQ(ierr);

  ierr = VecDuplicateVecs(svd->U[0],svd->n,&U);CHKERRQ(ierr);

  ierr = VecCopy(svd->vec_initial,V[0]);CHKERRQ(ierr);

  ierr = VecNormalize(V[0],&norm1);CHKERRQ(ierr);

  while (svd->reason == SVD_CONVERGED_ITERATING) {

    svd->its++;

    for (i=svd->nconv;i<svd->n;i++) {

      ierr = MatMult(svd->A,V[i],U[i]);CHKERRQ(ierr);

//      if (i>svd->nconv) { ierr = VecAXPY(U[i],-beta[i-1],U[i-1]);CHKERRQ(ierr); }

      ierr = cgs(U[i],i,U);CHKERRQ(ierr);

      ierr = cgs(U[i],i,U);CHKERRQ(ierr);

      ierr = VecNormalize(U[i],alpha+i-svd->nconv);CHKERRQ(ierr);

      if (svd->AT) {

        ierr = MatMult(svd->AT,U[i],V[i+1]);CHKERRQ(ierr);

      } else {

        ierr = MatMultTranspose(svd->A,U[i],V[i+1]);CHKERRQ(ierr);

      }

//      ierr = VecAXPY(V[i+1],-alpha[i],V[i]);CHKERRQ(ierr);

      ierr = cgs(V[i+1],i+1,V);CHKERRQ(ierr);

      ierr = cgs(V[i+1],i+1,V);CHKERRQ(ierr);

      ierr = VecNormalize(V[i+1],beta+i-svd->nconv);CHKERRQ(ierr);    

    }

    n = svd->n - svd->nconv;

    ierr = PetscMemzero(PT,sizeof(PetscScalar)*n*n);CHKERRQ(ierr);

    ierr = PetscMemzero(Q,sizeof(PetscScalar)*n*n);CHKERRQ(ierr);

    for (i=0;i<n;i++)

      PT[i*n+i] = Q[i*n+i] = 1.0;

    LAPACKbdsqr_("U",&n,&n,&n,&zero,alpha,beta,PT,&n,Q,&n,PETSC_NULL,&n,work,&info,1);

    k = svd->nconv;

    for (i=svd->nconv;i<svd->n;i++) {

      if (svd->which == SVD_SMALLEST) j = n-i+svd->nconv-1;

      else j = i-svd->nconv;

      svd->sigma[i] = alpha[j];

      ierr = VecSet(svd->V[i],0.0);CHKERRQ(ierr);

      for (l=0;l<n;l++) {

        ierr = VecAXPY(svd->V[i],PT[l*n+j],V[l+svd->nconv]);CHKERRQ(ierr);

      }

      ierr = VecSet(svd->U[i],0.0);CHKERRQ(ierr);

      ierr = VecMAXPY(svd->U[i],n,Q+j*n,U+svd->nconv);CHKERRQ(ierr);

      ierr = computeres(svd,svd->sigma[i],svd->U[i],svd->V[i],&norm1,&norm2);CHKERRQ(ierr);

      svd->errest[i] = sqrt(norm1*norm1+norm2*norm2);

      if (svd->errest[i] < svd->tol) {

        k++;

      } else break;

    }

    for (i=svd->nconv;i<k;i++) {

      ierr = VecCopy(svd->V[i],V[i]);CHKERRQ(ierr);

      ierr = VecCopy(svd->U[i],U[i]);CHKERRQ(ierr);

    }

    svd->nconv = k;

    SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,svd->ncv);

    if (svd->its > svd->max_it) svd->reason = SVD_DIVERGED_ITS;

    if (svd->nconv >= svd->nsv) svd->reason = SVD_CONVERGED_TOL;

    if (svd->reason == SVD_CONVERGED_ITERATING) {

      ierr = VecCopy(svd->V[svd->nconv],V[svd->nconv]);CHKERRQ(ierr);

    }

  }

  ierr = VecDestroyVecs(V,n+1);CHKERRQ(ierr);

  ierr = VecDestroyVecs(U,n);CHKERRQ(ierr);

  ierr = PetscFree(alpha);CHKERRQ(ierr);

  ierr = PetscFree(beta);CHKERRQ(ierr);

  ierr = PetscFree(Q);CHKERRQ(ierr);

  ierr = PetscFree(PT);CHKERRQ(ierr);

  ierr = PetscFree(work);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

EXTERN_C_BEGIN

#undef __FUNCT__  

#define __FUNCT__ "SVDCreate_LANCZOS"

PetscErrorCode SVDCreate_LANCZOS(SVD svd)

{

  PetscFunctionBegin;

  svd->ops->setup = SVDSetUp_LANCZOS;

  svd->ops->solve = SVDSolve_LANCZOS;

  PetscFunctionReturn(0);

}

EXTERN_C_END