Actual source code: svdelemen.cxx

slepc-3.23.0 2025-03-29
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    This file implements a wrapper to the Elemental SVD solver
 12: */

 14: #include <slepc/private/svdimpl.h>
 15: #include <petsc/private/petscelemental.h>

 17: typedef struct {
 18:   Mat Ae;        /* converted matrix */
 19: } SVD_Elemental;

 21: static PetscErrorCode SVDSetUp_Elemental(SVD svd)
 22: {
 23:   SVD_Elemental  *ctx = (SVD_Elemental*)svd->data;
 24:   PetscInt       M,N;

 26:   PetscFunctionBegin;
 27:   SVDCheckStandard(svd);
 28:   SVDCheckDefinite(svd);
 29:   if (svd->nsv==0) svd->nsv = 1;
 30:   PetscCall(MatGetSize(svd->A,&M,&N));
 31:   PetscCheck(M==N,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Not implemented for rectangular matrices");
 32:   svd->ncv = N;
 33:   if (svd->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(svd,"Warning: parameter mpd ignored\n"));
 34:   if (svd->max_it==PETSC_DETERMINE) svd->max_it = 1;
 35:   svd->leftbasis = PETSC_TRUE;
 36:   SVDCheckIgnored(svd,SVD_FEATURE_STOPPING);
 37:   PetscCall(SVDAllocateSolution(svd,0));

 39:   /* convert matrix */
 40:   PetscCall(MatDestroy(&ctx->Ae));
 41:   PetscCall(MatConvert(svd->OP,MATELEMENTAL,MAT_INITIAL_MATRIX,&ctx->Ae));
 42:   PetscFunctionReturn(PETSC_SUCCESS);
 43: }

 45: static PetscErrorCode SVDSolve_Elemental(SVD svd)
 46: {
 47:   SVD_Elemental  *ctx = (SVD_Elemental*)svd->data;
 48:   Mat            A = ctx->Ae,Z,Q,U,V;
 49:   Mat_Elemental  *a = (Mat_Elemental*)A->data,*q,*z;
 50:   PetscInt       i,rrank,ridx,erow;

 52:   PetscFunctionBegin;
 53:   El::DistMatrix<PetscReal,El::STAR,El::VC> sigma(*a->grid);
 54:   PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Z));
 55:   PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Q));
 56:   z = (Mat_Elemental*)Z->data;
 57:   q = (Mat_Elemental*)Q->data;

 59:   El::SVD(*a->emat,*z->emat,sigma,*q->emat);

 61:   for (i=0;i<svd->ncv;i++) {
 62:     P2RO(A,1,i,&rrank,&ridx);
 63:     RO2E(A,1,rrank,ridx,&erow);
 64:     svd->sigma[i] = sigma.Get(erow,0);
 65:   }
 66:   PetscCall(BVGetMat(svd->U,&U));
 67:   PetscCall(MatConvert(Z,MATDENSE,MAT_REUSE_MATRIX,&U));
 68:   PetscCall(BVRestoreMat(svd->U,&U));
 69:   PetscCall(MatDestroy(&Z));
 70:   PetscCall(BVGetMat(svd->V,&V));
 71:   PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&V));
 72:   PetscCall(BVRestoreMat(svd->V,&V));
 73:   PetscCall(MatDestroy(&Q));

 75:   svd->nconv  = svd->ncv;
 76:   svd->its    = 1;
 77:   svd->reason = SVD_CONVERGED_TOL;
 78:   PetscFunctionReturn(PETSC_SUCCESS);
 79: }

 81: static PetscErrorCode SVDDestroy_Elemental(SVD svd)
 82: {
 83:   PetscFunctionBegin;
 84:   PetscCall(PetscFree(svd->data));
 85:   PetscFunctionReturn(PETSC_SUCCESS);
 86: }

 88: static PetscErrorCode SVDReset_Elemental(SVD svd)
 89: {
 90:   SVD_Elemental  *ctx = (SVD_Elemental*)svd->data;

 92:   PetscFunctionBegin;
 93:   PetscCall(MatDestroy(&ctx->Ae));
 94:   PetscFunctionReturn(PETSC_SUCCESS);
 95: }

 97: SLEPC_EXTERN PetscErrorCode SVDCreate_Elemental(SVD svd)
 98: {
 99:   SVD_Elemental  *ctx;

101:   PetscFunctionBegin;
102:   PetscCall(PetscNew(&ctx));
103:   svd->data = (void*)ctx;

105:   svd->ops->solve          = SVDSolve_Elemental;
106:   svd->ops->setup          = SVDSetUp_Elemental;
107:   svd->ops->destroy        = SVDDestroy_Elemental;
108:   svd->ops->reset          = SVDReset_Elemental;
109:   PetscFunctionReturn(PETSC_SUCCESS);
110: }