Subversion Repositories slepc-dev

/*                      

/*                      

   SLEPc singular value solver: "cross"

   SLEPc singular value solver: "cross"

   Method: Uses a Hermitian eigensolver for A^T*A

   Method: Uses a Hermitian eigensolver for A^T*A

   Last update: Jan 2007

   Last update: Jan 2007

*/

*/

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

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

#include "slepceps.h"

#include "slepceps.h"

typedef struct {

typedef struct {

  EPS eps;

  EPS eps;

  Mat mat;

  Mat mat;

  Vec w,diag;

  Vec w,diag;

} SVD_CROSS;

} SVD_CROSS;

#undef __FUNCT__  

#undef __FUNCT__  

#define __FUNCT__ "ShellMatMult_CROSS"

#define __FUNCT__ "ShellMatMult_CROSS"

PetscErrorCode ShellMatMult_CROSS(Mat B,Vec x, Vec y)

PetscErrorCode ShellMatMult_CROSS(Mat B,Vec x, Vec y)

{

{

  PetscErrorCode ierr;

  PetscErrorCode ierr;

  SVD            svd;

  SVD            svd;

  SVD_CROSS      *cross;

  SVD_CROSS      *cross;

  PetscFunctionBegin;

  PetscFunctionBegin;

  ierr = MatShellGetContext(B,(void**)&svd);CHKERRQ(ierr);

  ierr = MatShellGetContext(B,(void**)&svd);CHKERRQ(ierr);

  cross = (SVD_CROSS *)svd->data;

  cross = (SVD_CROSS *)svd->data;

  ierr = SVDMatMult(svd,PETSC_FALSE,x,cross->w);CHKERRQ(ierr);

  ierr = SVDMatMult(svd,PETSC_FALSE,x,cross->w);CHKERRQ(ierr);

  ierr = SVDMatMult(svd,PETSC_TRUE,cross->w,y);CHKERRQ(ierr);

  ierr = SVDMatMult(svd,PETSC_TRUE,cross->w,y);CHKERRQ(ierr);

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

}

}

#undef __FUNCT__  

#undef __FUNCT__  

#define __FUNCT__ "ShellMatGetDiagonal_CROSS"

#define __FUNCT__ "ShellMatGetDiagonal_CROSS"

PetscErrorCode ShellMatGetDiagonal_CROSS(Mat B,Vec d)

PetscErrorCode ShellMatGetDiagonal_CROSS(Mat B,Vec d)

{

{

  PetscErrorCode    ierr;

  PetscErrorCode    ierr;

  SVD               svd;

  SVD               svd;

  SVD_CROSS         *cross;

  SVD_CROSS         *cross;

  PetscInt          N,n,i,j,start,end,ncols;

  PetscInt          N,n,i,j,start,end,ncols;

  PetscScalar       *work1,*work2,*diag;

  PetscScalar       *work1,*work2,*diag;

  const PetscInt    *cols;

  const PetscInt    *cols;

  const PetscScalar *vals;

  const PetscScalar *vals;

  PetscFunctionBegin;

  PetscFunctionBegin;

  ierr = MatShellGetContext(B,(void**)&svd);CHKERRQ(ierr);

  ierr = MatShellGetContext(B,(void**)&svd);CHKERRQ(ierr);

  cross = (SVD_CROSS *)svd->data;

  cross = (SVD_CROSS *)svd->data;

  if (!cross->diag) {

  if (!cross->diag) {

    /* compute diagonal from rows and store in cross->diag */

    /* compute diagonal from rows and store in cross->diag */

    ierr = VecDuplicate(d,&cross->diag);CHKERRQ(ierr);

    ierr = VecDuplicate(d,&cross->diag);CHKERRQ(ierr);

    ierr = SVDMatGetSize(svd,PETSC_NULL,&N);CHKERRQ(ierr);

    ierr = SVDMatGetSize(svd,PETSC_NULL,&N);CHKERRQ(ierr);

    ierr = SVDMatGetLocalSize(svd,PETSC_NULL,&n);CHKERRQ(ierr);

    ierr = SVDMatGetLocalSize(svd,PETSC_NULL,&n);CHKERRQ(ierr);

    ierr = PetscMalloc(sizeof(PetscScalar)*N,&work1);CHKERRQ(ierr);

    ierr = PetscMalloc(sizeof(PetscScalar)*N,&work1);CHKERRQ(ierr);

    ierr = PetscMalloc(sizeof(PetscScalar)*N,&work2);CHKERRQ(ierr);

    ierr = PetscMalloc(sizeof(PetscScalar)*N,&work2);CHKERRQ(ierr);

    for (i=0;i<n;i++) work1[i] = work2[i] = 0.0;

    for (i=0;i<n;i++) work1[i] = work2[i] = 0.0;

    if (svd->AT) {

    if (svd->AT) {

      ierr = MatGetOwnershipRange(svd->AT,&start,&end);CHKERRQ(ierr);

      ierr = MatGetOwnershipRange(svd->AT,&start,&end);CHKERRQ(ierr);

      for (i=start;i<end;i++) {

      for (i=start;i<end;i++) {

        ierr = MatGetRow(svd->AT,i,&ncols,PETSC_NULL,&vals);CHKERRQ(ierr);

        ierr = MatGetRow(svd->AT,i,&ncols,PETSC_NULL,&vals);CHKERRQ(ierr);

        for (j=0;j<ncols;j++)

        for (j=0;j<ncols;j++)

          work1[i] += vals[j]*vals[j];

          work1[i] += vals[j]*vals[j];

        ierr = MatRestoreRow(svd->AT,i,&ncols,PETSC_NULL,&vals);CHKERRQ(ierr);

        ierr = MatRestoreRow(svd->AT,i,&ncols,PETSC_NULL,&vals);CHKERRQ(ierr);

      }

      }

    } else {

    } else {

      ierr = MatGetOwnershipRange(svd->A,&start,&end);CHKERRQ(ierr);

      ierr = MatGetOwnershipRange(svd->A,&start,&end);CHKERRQ(ierr);

      for (i=start;i<end;i++) {

      for (i=start;i<end;i++) {

        ierr = MatGetRow(svd->A,i,&ncols,&cols,&vals);CHKERRQ(ierr);

        ierr = MatGetRow(svd->A,i,&ncols,&cols,&vals);CHKERRQ(ierr);

        for (j=0;j<ncols;j++)

        for (j=0;j<ncols;j++)

          work1[cols[j]] += vals[j]*vals[j];

          work1[cols[j]] += vals[j]*vals[j];

        ierr = MatRestoreRow(svd->A,i,&ncols,&cols,&vals);CHKERRQ(ierr);

        ierr = MatRestoreRow(svd->A,i,&ncols,&cols,&vals);CHKERRQ(ierr);

      }

      }

    }

    }

    ierr = MPI_Allreduce(work1,work2,N,MPIU_SCALAR,MPI_SUM,svd->comm);CHKERRQ(ierr);

    ierr = MPI_Allreduce(work1,work2,N,MPIU_SCALAR,MPI_SUM,svd->comm);CHKERRQ(ierr);

    ierr = VecGetOwnershipRange(cross->diag,&start,&end);CHKERRQ(ierr);

    ierr = VecGetOwnershipRange(cross->diag,&start,&end);CHKERRQ(ierr);

    ierr = VecGetArray(cross->diag,&diag);CHKERRQ(ierr);

    ierr = VecGetArray(cross->diag,&diag);CHKERRQ(ierr);

    for (i=start;i<end;i++)

    for (i=start;i<end;i++)

      diag[i-start] = work2[i];

      diag[i-start] = work2[i];

    ierr = VecRestoreArray(cross->diag,&diag);CHKERRQ(ierr);

    ierr = VecRestoreArray(cross->diag,&diag);CHKERRQ(ierr);

    ierr = PetscFree(work1);CHKERRQ(ierr);

    ierr = PetscFree(work1);CHKERRQ(ierr);

    ierr = PetscFree(work2);CHKERRQ(ierr);

    ierr = PetscFree(work2);CHKERRQ(ierr);

  }

  }

  ierr = VecCopy(cross->diag,d);CHKERRQ(ierr);

  ierr = VecCopy(cross->diag,d);CHKERRQ(ierr);

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

}

}

#undef __FUNCT__  

#undef __FUNCT__  

#define __FUNCT__ "SVDSetUp_CROSS"

#define __FUNCT__ "SVDSetUp_CROSS"

PetscErrorCode SVDSetUp_CROSS(SVD svd)

PetscErrorCode SVDSetUp_CROSS(SVD svd)

{

{

  PetscErrorCode    ierr;

  PetscErrorCode    ierr;

  SVD_CROSS         *cross = (SVD_CROSS *)svd->data;

  SVD_CROSS         *cross = (SVD_CROSS *)svd->data;

  PetscInt          n;

  PetscInt          n;

  PetscFunctionBegin;

  PetscFunctionBegin;

  if (cross->mat) {

  if (cross->mat) {

     ierr = MatDestroy(cross->mat);CHKERRQ(ierr);

     ierr = MatDestroy(cross->mat);CHKERRQ(ierr);

     ierr = VecDestroy(cross->w);CHKERRQ(ierr);

     ierr = VecDestroy(cross->w);CHKERRQ(ierr);

  }

  }

  if (cross->diag) {

  if (cross->diag) {

     ierr = VecDestroy(cross->diag);CHKERRQ(ierr);

     ierr = VecDestroy(cross->diag);CHKERRQ(ierr);

  }

  }

  ierr = SVDMatGetLocalSize(svd,PETSC_NULL,&n);CHKERRQ(ierr);

  ierr = SVDMatGetLocalSize(svd,PETSC_NULL,&n);CHKERRQ(ierr);

  ierr = MatCreateShell(svd->comm,n,n,PETSC_DETERMINE,PETSC_DETERMINE,svd,&cross->mat);CHKERRQ(ierr);

  ierr = MatCreateShell(svd->comm,n,n,PETSC_DETERMINE,PETSC_DETERMINE,svd,&cross->mat);CHKERRQ(ierr);

  ierr = MatShellSetOperation(cross->mat,MATOP_MULT,(void(*)(void))ShellMatMult_CROSS);CHKERRQ(ierr);  

  ierr = MatShellSetOperation(cross->mat,MATOP_MULT,(void(*)(void))ShellMatMult_CROSS);CHKERRQ(ierr);  

  ierr = MatShellSetOperation(cross->mat,MATOP_GET_DIAGONAL,(void(*)(void))ShellMatGetDiagonal_CROSS);CHKERRQ(ierr);  

  ierr = MatShellSetOperation(cross->mat,MATOP_GET_DIAGONAL,(void(*)(void))ShellMatGetDiagonal_CROSS);CHKERRQ(ierr);  

  ierr = SVDMatGetVecs(svd,PETSC_NULL,&cross->w);CHKERRQ(ierr);

  ierr = SVDMatGetVecs(svd,PETSC_NULL,&cross->w);CHKERRQ(ierr);

  ierr = EPSSetOperators(cross->eps,cross->mat,PETSC_NULL);CHKERRQ(ierr);

  ierr = EPSSetOperators(cross->eps,cross->mat,PETSC_NULL);CHKERRQ(ierr);

  ierr = EPSSetProblemType(cross->eps,EPS_HEP);CHKERRQ(ierr);

  ierr = EPSSetProblemType(cross->eps,EPS_HEP);CHKERRQ(ierr);

  ierr = EPSSetDimensions(cross->eps,svd->nsv,svd->ncv);CHKERRQ(ierr);

  ierr = EPSSetDimensions(cross->eps,svd->nsv,svd->ncv);CHKERRQ(ierr);

  ierr = EPSSetTolerances(cross->eps,svd->tol,svd->max_it);CHKERRQ(ierr);

  ierr = EPSSetTolerances(cross->eps,svd->tol,svd->max_it);CHKERRQ(ierr);

  ierr = EPSSetUp(cross->eps);CHKERRQ(ierr);

  ierr = EPSSetUp(cross->eps);CHKERRQ(ierr);

  ierr = EPSGetDimensions(cross->eps,PETSC_NULL,&svd->ncv);CHKERRQ(ierr);

  ierr = EPSGetDimensions(cross->eps,PETSC_NULL,&svd->ncv);CHKERRQ(ierr);

  ierr = EPSGetTolerances(cross->eps,&svd->tol,&svd->max_it);CHKERRQ(ierr);

  ierr = EPSGetTolerances(cross->eps,&svd->tol,&svd->max_it);CHKERRQ(ierr);

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

}

}

#undef __FUNCT__  

#undef __FUNCT__  

#define __FUNCT__ "SVDSolve_CROSS"

#define __FUNCT__ "SVDSolve_CROSS"

PetscErrorCode SVDSolve_CROSS(SVD svd)

PetscErrorCode SVDSolve_CROSS(SVD svd)

{

{

  PetscErrorCode ierr;

  PetscErrorCode ierr;

  SVD_CROSS      *cross = (SVD_CROSS *)svd->data;

  SVD_CROSS      *cross = (SVD_CROSS *)svd->data;

  int            i;

  int            i;

  PetscScalar    sigma;

  PetscScalar    sigma;

  PetscFunctionBegin;

  PetscFunctionBegin;

  ierr = EPSSetWhichEigenpairs(cross->eps,svd->which == SVD_LARGEST ? EPS_LARGEST_REAL : EPS_SMALLEST_MAGNITUDE);CHKERRQ(ierr);

  ierr = EPSSetWhichEigenpairs(cross->eps,svd->which == SVD_LARGEST ? EPS_LARGEST_REAL : EPS_SMALLEST_MAGNITUDE);CHKERRQ(ierr);

  ierr = EPSSetInitialVector(cross->eps,svd->vec_initial);CHKERRQ(ierr);

  ierr = EPSSetInitialVector(cross->eps,svd->vec_initial);CHKERRQ(ierr);

  ierr = EPSSolve(cross->eps);CHKERRQ(ierr);

  ierr = EPSSolve(cross->eps);CHKERRQ(ierr);

  ierr = EPSGetConverged(cross->eps,&svd->nconv);CHKERRQ(ierr);

  ierr = EPSGetConverged(cross->eps,&svd->nconv);CHKERRQ(ierr);

  ierr = EPSGetIterationNumber(cross->eps,&svd->its);CHKERRQ(ierr);

  ierr = EPSGetIterationNumber(cross->eps,&svd->its);CHKERRQ(ierr);

  ierr = EPSGetConvergedReason(cross->eps,(EPSConvergedReason*)&svd->reason);CHKERRQ(ierr);

  ierr = EPSGetConvergedReason(cross->eps,(EPSConvergedReason*)&svd->reason);CHKERRQ(ierr);

  ierr = EPSGetOperationCounters(cross->eps,PETSC_NULL,&svd->dots,PETSC_NULL);CHKERRQ(ierr);

  ierr = EPSGetOperationCounters(cross->eps,PETSC_NULL,&svd->dots,PETSC_NULL);CHKERRQ(ierr);

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

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

    ierr = EPSGetEigenpair(cross->eps,i,&sigma,PETSC_NULL,svd->V[i],PETSC_NULL);CHKERRQ(ierr);

    ierr = EPSGetEigenpair(cross->eps,i,&sigma,PETSC_NULL,svd->V[i],PETSC_NULL);CHKERRQ(ierr);

    svd->sigma[i] = sqrt(PetscRealPart(sigma));

    svd->sigma[i] = sqrt(PetscRealPart(sigma));

  }

  }

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

}

}

#undef __FUNCT__  

#undef __FUNCT__  

#define __FUNCT__ "SVDMonitor_CROSS"

#define __FUNCT__ "SVDMonitor_CROSS"

PetscErrorCode SVDMonitor_CROSS(EPS eps,int its,int nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,int nest,void *ctx)

PetscErrorCode SVDMonitor_CROSS(EPS eps,int its,int nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,int nest,void *ctx)

{

{

  int       i;

  int       i;

  SVD       svd = (SVD)ctx;

  SVD       svd = (SVD)ctx;

  PetscFunctionBegin;

  PetscFunctionBegin;

  for (i=0;i<nest;i++) {

  for (i=0;i<nest;i++) {

    svd->sigma[i] = sqrt(PetscRealPart(eigr[i]));

    svd->sigma[i] = sqrt(PetscRealPart(eigr[i]));

    svd->errest[i] = errest[i];

    svd->errest[i] = errest[i];

  }

  }

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

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

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

}

}

#undef __FUNCT__  

#undef __FUNCT__  

#define __FUNCT__ "SVDSetFromOptions_CROSS"

#define __FUNCT__ "SVDSetFromOptions_CROSS"

PetscErrorCode SVDSetFromOptions_CROSS(SVD svd)

PetscErrorCode SVDSetFromOptions_CROSS(SVD svd)

{

{

  PetscErrorCode ierr;

  PetscErrorCode ierr;

  SVD_CROSS      *cross = (SVD_CROSS *)svd->data;

  SVD_CROSS      *cross = (SVD_CROSS *)svd->data;

  PetscFunctionBegin;

  PetscFunctionBegin;

  ierr = EPSSetFromOptions(cross->eps);CHKERRQ(ierr);

  ierr = EPSSetFromOptions(cross->eps);CHKERRQ(ierr);

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

}

}

EXTERN_C_BEGIN

EXTERN_C_BEGIN

#undef __FUNCT__  

#undef __FUNCT__  

#define __FUNCT__ "SVDCrossSetEPS_CROSS"

#define __FUNCT__ "SVDCrossSetEPS_CROSS"

PetscErrorCode SVDCrossSetEPS_CROSS(SVD svd,EPS eps)

PetscErrorCode SVDCrossSetEPS_CROSS(SVD svd,EPS eps)

{

{

  PetscErrorCode ierr;

  PetscErrorCode ierr;

  SVD_CROSS      *cross = (SVD_CROSS *)svd->data;

  SVD_CROSS      *cross = (SVD_CROSS *)svd->data;

  PetscFunctionBegin;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE,2);

  PetscValidHeaderSpecific(eps,EPS_COOKIE,2);

  PetscCheckSameComm(svd,1,eps,2);

  PetscCheckSameComm(svd,1,eps,2);

  ierr = PetscObjectReference((PetscObject)eps);CHKERRQ(ierr);

  ierr = PetscObjectReference((PetscObject)eps);CHKERRQ(ierr);

  ierr = EPSDestroy(cross->eps);CHKERRQ(ierr);  

  ierr = EPSDestroy(cross->eps);CHKERRQ(ierr);  

  cross->eps = eps;

  cross->eps = eps;

  svd->setupcalled = 0;

  svd->setupcalled = 0;

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

}

}

EXTERN_C_END

EXTERN_C_END

#undef __FUNCT__  

#undef __FUNCT__  

#define __FUNCT__ "SVDCrossSetEPS"

#define __FUNCT__ "SVDCrossSetEPS"

/*@

/*@

   SVDCrossSetEPS - Associate an eigensolver object (EPS) to the

   SVDCrossSetEPS - Associate an eigensolver object (EPS) to the

   singular value solver.

   singular value solver.

   Collective on SVD

   Collective on SVD

   Input Parameters:

   Input Parameters:

+  svd - singular value solver

+  svd - singular value solver

-  eps - the eigensolver object

-  eps - the eigensolver object

   Level: advanced

   Level: advanced

.seealso: SVDCrossGetEPS()

.seealso: SVDCrossGetEPS()

@*/

@*/

PetscErrorCode SVDCrossSetEPS(SVD svd,EPS eps)

PetscErrorCode SVDCrossSetEPS(SVD svd,EPS eps)

{

{

  PetscErrorCode ierr, (*f)(SVD,EPS eps);

  PetscErrorCode ierr, (*f)(SVD,EPS eps);

  PetscFunctionBegin;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(svd,SVD_COOKIE,1);

  PetscValidHeaderSpecific(svd,SVD_COOKIE,1);

  ierr = PetscObjectQueryFunction((PetscObject)svd,"SVDCrossSetEPS_C",(void (**)())&f);CHKERRQ(ierr);

  ierr = PetscObjectQueryFunction((PetscObject)svd,"SVDCrossSetEPS_C",(void (**)())&f);CHKERRQ(ierr);

  if (f) {

  if (f) {

    ierr = (*f)(svd,eps);CHKERRQ(ierr);

    ierr = (*f)(svd,eps);CHKERRQ(ierr);

  }

  }

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

}

}

EXTERN_C_BEGIN

EXTERN_C_BEGIN

#undef __FUNCT__  

#undef __FUNCT__  

#define __FUNCT__ "SVDCrossGetEPS_CROSS"

#define __FUNCT__ "SVDCrossGetEPS_CROSS"

PetscErrorCode SVDCrossGetEPS_CROSS(SVD svd,EPS *eps)

PetscErrorCode SVDCrossGetEPS_CROSS(SVD svd,EPS *eps)

{

{

  SVD_CROSS *cross = (SVD_CROSS *)svd->data;

  SVD_CROSS *cross = (SVD_CROSS *)svd->data;

  PetscFunctionBegin;

  PetscFunctionBegin;

  PetscValidPointer(eps,2);

  PetscValidPointer(eps,2);

  *eps = cross->eps;

  *eps = cross->eps;

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

}

}

EXTERN_C_END

EXTERN_C_END

#undef __FUNCT__  

#undef __FUNCT__  

#define __FUNCT__ "SVDCrossGetEPS"

#define __FUNCT__ "SVDCrossGetEPS"

/*@C

/*@C

   SVDCrossGetEPS - Retrieve the eigensolver object (EPS) associated

   SVDCrossGetEPS - Retrieve the eigensolver object (EPS) associated

   to the singular value solver.

   to the singular value solver.

   Not Collective

   Not Collective

   Input Parameter:

   Input Parameter:

.  svd - singular value solver

.  svd - singular value solver

   Output Parameter:

   Output Parameter:

.  eps - the eigensolver object

.  eps - the eigensolver object

   Level: advanced

   Level: advanced

.seealso: SVDCrossSetEPS()

.seealso: SVDCrossSetEPS()

@*/

@*/

PetscErrorCode SVDCrossGetEPS(SVD svd,EPS *eps)

PetscErrorCode SVDCrossGetEPS(SVD svd,EPS *eps)

{

{

  PetscErrorCode ierr, (*f)(SVD,EPS *eps);

  PetscErrorCode ierr, (*f)(SVD,EPS *eps);

  PetscFunctionBegin;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(svd,SVD_COOKIE,1);

  PetscValidHeaderSpecific(svd,SVD_COOKIE,1);

  ierr = PetscObjectQueryFunction((PetscObject)svd,"SVDCrossGetEPS_C",(void (**)())&f);CHKERRQ(ierr);

  ierr = PetscObjectQueryFunction((PetscObject)svd,"SVDCrossGetEPS_C",(void (**)())&f);CHKERRQ(ierr);

  if (f) {

  if (f) {

    ierr = (*f)(svd,eps);CHKERRQ(ierr);

    ierr = (*f)(svd,eps);CHKERRQ(ierr);

  }

  }

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

}

}

#undef __FUNCT__  

#undef __FUNCT__  

#define __FUNCT__ "SVDView_CROSS"

#define __FUNCT__ "SVDView_CROSS"

PetscErrorCode SVDView_CROSS(SVD svd,PetscViewer viewer)

PetscErrorCode SVDView_CROSS(SVD svd,PetscViewer viewer)

{

{

  PetscErrorCode ierr;

  PetscErrorCode ierr;

  SVD_CROSS      *cross = (SVD_CROSS *)svd->data;

  SVD_CROSS      *cross = (SVD_CROSS *)svd->data;

  PetscFunctionBegin;

  PetscFunctionBegin;

  ierr = EPSView(cross->eps,viewer);CHKERRQ(ierr);

  ierr = EPSView(cross->eps,viewer);CHKERRQ(ierr);

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

}

}

#undef __FUNCT__  

#undef __FUNCT__  

#define __FUNCT__ "SVDDestroy_CROSS"

#define __FUNCT__ "SVDDestroy_CROSS"

PetscErrorCode SVDDestroy_CROSS(SVD svd)

PetscErrorCode SVDDestroy_CROSS(SVD svd)

{

{

  PetscErrorCode ierr;

  PetscErrorCode ierr;

  SVD_CROSS      *cross = (SVD_CROSS *)svd->data;

  SVD_CROSS      *cross = (SVD_CROSS *)svd->data;

  PetscFunctionBegin;

  PetscFunctionBegin;

  ierr = EPSDestroy(cross->eps);CHKERRQ(ierr);

  ierr = EPSDestroy(cross->eps);CHKERRQ(ierr);

  if (cross->mat) {

  if (cross->mat) {

    ierr = MatDestroy(cross->mat);CHKERRQ(ierr);

    ierr = MatDestroy(cross->mat);CHKERRQ(ierr);

    ierr = VecDestroy(cross->w);CHKERRQ(ierr);

    ierr = VecDestroy(cross->w);CHKERRQ(ierr);

  }

  }

  if (cross->diag) {

  if (cross->diag) {

    ierr = VecDestroy(cross->diag);CHKERRQ(ierr);

    ierr = VecDestroy(cross->diag);CHKERRQ(ierr);

  }

  }

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

}

}

EXTERN_C_BEGIN

EXTERN_C_BEGIN

#undef __FUNCT__  

#undef __FUNCT__  

#define __FUNCT__ "SVDCreate_CROSS"

#define __FUNCT__ "SVDCreate_CROSS"

PetscErrorCode SVDCreate_CROSS(SVD svd)

PetscErrorCode SVDCreate_CROSS(SVD svd)

{

{

  PetscErrorCode ierr;

  PetscErrorCode ierr;

  SVD_CROSS      *cross;

  SVD_CROSS      *cross;

  ST             st;

  ST             st;

  PetscFunctionBegin;

  PetscFunctionBegin;

  ierr = PetscNew(SVD_CROSS,&cross);CHKERRQ(ierr);

  ierr = PetscNew(SVD_CROSS,&cross);CHKERRQ(ierr);

  PetscLogObjectMemory(svd,sizeof(SVD_CROSS));

  PetscLogObjectMemory(svd,sizeof(SVD_CROSS));

  svd->data                = (void *)cross;

  svd->data                = (void *)cross;

  svd->ops->solve          = SVDSolve_CROSS;

  svd->ops->solve          = SVDSolve_CROSS;

  svd->ops->setup          = SVDSetUp_CROSS;

  svd->ops->setup          = SVDSetUp_CROSS;

  svd->ops->setfromoptions = SVDSetFromOptions_CROSS;

  svd->ops->setfromoptions = SVDSetFromOptions_CROSS;

  svd->ops->destroy        = SVDDestroy_CROSS;

  svd->ops->destroy        = SVDDestroy_CROSS;

  svd->ops->view           = SVDView_CROSS;

  svd->ops->view           = SVDView_CROSS;

  ierr = PetscObjectComposeFunctionDynamic((PetscObject)svd,"SVDCrossSetEPS_C","SVDCrossSetEPS_CROSS",SVDCrossSetEPS_CROSS);CHKERRQ(ierr);

  ierr = PetscObjectComposeFunctionDynamic((PetscObject)svd,"SVDCrossSetEPS_C","SVDCrossSetEPS_CROSS",SVDCrossSetEPS_CROSS);CHKERRQ(ierr);

  ierr = PetscObjectComposeFunctionDynamic((PetscObject)svd,"SVDCrossGetEPS_C","SVDCrossGetEPS_CROSS",SVDCrossGetEPS_CROSS);CHKERRQ(ierr);

  ierr = PetscObjectComposeFunctionDynamic((PetscObject)svd,"SVDCrossGetEPS_C","SVDCrossGetEPS_CROSS",SVDCrossGetEPS_CROSS);CHKERRQ(ierr);

  ierr = EPSCreate(svd->comm,&cross->eps);CHKERRQ(ierr);

  ierr = EPSCreate(svd->comm,&cross->eps);CHKERRQ(ierr);

  ierr = EPSSetOptionsPrefix(cross->eps,svd->prefix);CHKERRQ(ierr);

  ierr = EPSSetOptionsPrefix(cross->eps,svd->prefix);CHKERRQ(ierr);

  ierr = EPSAppendOptionsPrefix(cross->eps,"svd_");CHKERRQ(ierr);

  ierr = EPSAppendOptionsPrefix(cross->eps,"svd_");CHKERRQ(ierr);

  PetscLogObjectParent(svd,cross->eps);

  PetscLogObjectParent(svd,cross->eps);

  ierr = EPSSetWhichEigenpairs(cross->eps,EPS_LARGEST_REAL);CHKERRQ(ierr);

  ierr = EPSSetWhichEigenpairs(cross->eps,EPS_LARGEST_REAL);CHKERRQ(ierr);

  ierr = EPSMonitorSet(cross->eps,SVDMonitor_CROSS,svd,PETSC_NULL);CHKERRQ(ierr);

  ierr = EPSMonitorSet(cross->eps,SVDMonitor_CROSS,svd,PETSC_NULL);CHKERRQ(ierr);

  ierr = EPSGetST(cross->eps,&st);CHKERRQ(ierr);

  ierr = EPSGetST(cross->eps,&st);CHKERRQ(ierr);

  ierr = STSetMatMode(st,STMATMODE_SHELL);CHKERRQ(ierr);

  ierr = STSetMatMode(st,STMATMODE_SHELL);CHKERRQ(ierr);

  cross->mat = PETSC_NULL;

  cross->mat = PETSC_NULL;

  cross->w = PETSC_NULL;

  cross->w = PETSC_NULL;

  cross->diag = PETSC_NULL;

  cross->diag = PETSC_NULL;

  PetscFunctionReturn(0);

  PetscFunctionReturn(0);

}

}

EXTERN_C_END

EXTERN_C_END