Actual source code: svdelemen.cxx
slepc-3.23.0 2025-03-29
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: }