Actual source code: elpa.c
slepc-3.23.1 2025-05-01
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 eigensolvers in ELPA.
12: */
14: #include <slepc/private/epsimpl.h>
15: #include <slepc/private/slepcscalapack.h>
16: #include <elpa/elpa.h>
18: #define PetscCallELPA(func, ...) do { \
19: PetscErrorCode elpa_ierr_; \
20: PetscStackPushExternal(PetscStringize(func)); \
21: func(__VA_ARGS__,&elpa_ierr_); \
22: PetscStackPop; \
23: PetscCheck(!elpa_ierr_,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error calling %s: error code %d",PetscStringize(func(__VA_ARGS__,&elpa_ierr)),elpa_ierr_); \
24: } while (0)
26: #define PetscCallELPARET(func, ...) do { \
27: PetscStackPushExternal(PetscStringize(func)); \
28: PetscErrorCode elpa_ierr_ = func(__VA_ARGS__); \
29: PetscStackPop; \
30: PetscCheck(!elpa_ierr_,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error calling %s: error code %d",PetscStringize(func(__VA_ARGS__)),elpa_ierr_); \
31: } while (0)
33: #define PetscCallELPANOARG(func) do { \
34: PetscErrorCode elpa_ierr_; \
35: PetscStackPushExternal(PetscStringize(func)); \
36: func(&elpa_ierr_); \
37: PetscStackPop; \
38: PetscCheck(!elpa_ierr_,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error calling %s: error code %d",PetscStringize(func(&elpa_ierr)),elpa_ierr_); \
39: } while (0)
41: typedef struct {
42: Mat As,Bs; /* converted matrices */
43: } EPS_ELPA;
45: static PetscErrorCode EPSSetUp_ELPA(EPS eps)
46: {
47: EPS_ELPA *ctx = (EPS_ELPA*)eps->data;
48: Mat A,B;
49: PetscInt nmat;
50: PetscBool isshift;
51: PetscScalar shift;
53: PetscFunctionBegin;
54: EPSCheckHermitianDefinite(eps);
55: EPSCheckNotStructured(eps);
56: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift));
57: PetscCheck(isshift,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support spectral transformations");
58: if (eps->nev==0) eps->nev = 1;
59: eps->ncv = eps->n;
60: if (eps->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
61: if (eps->max_it==PETSC_DETERMINE) eps->max_it = 1;
62: if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
63: PetscCheck(eps->which!=EPS_ALL || eps->inta==eps->intb,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support interval computation");
64: EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION);
65: EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION | EPS_FEATURE_CONVERGENCE | EPS_FEATURE_STOPPING);
66: PetscCall(EPSAllocateSolution(eps,0));
68: /* convert matrices */
69: PetscCall(MatDestroy(&ctx->As));
70: PetscCall(MatDestroy(&ctx->Bs));
71: PetscCall(STGetNumMatrices(eps->st,&nmat));
72: PetscCall(STGetMatrix(eps->st,0,&A));
73: PetscCall(MatConvert(A,MATSCALAPACK,MAT_INITIAL_MATRIX,&ctx->As));
74: if (nmat>1) {
75: PetscCall(STGetMatrix(eps->st,1,&B));
76: PetscCall(MatConvert(B,MATSCALAPACK,MAT_INITIAL_MATRIX,&ctx->Bs));
77: }
78: PetscCall(STGetShift(eps->st,&shift));
79: if (shift != 0.0) {
80: if (nmat>1) PetscCall(MatAXPY(ctx->As,-shift,ctx->Bs,SAME_NONZERO_PATTERN));
81: else PetscCall(MatShift(ctx->As,-shift));
82: }
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: static PetscErrorCode EPSSolve_ELPA(EPS eps)
87: {
88: EPS_ELPA *ctx = (EPS_ELPA*)eps->data;
89: Mat A = ctx->As,B = ctx->Bs,Q,V;
90: Mat_ScaLAPACK *a = (Mat_ScaLAPACK*)A->data,*b,*q;
91: PetscReal *w = eps->errest; /* used to store real eigenvalues */
92: PetscInt i;
93: elpa_t handle;
95: PetscFunctionBegin;
96: PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Q));
97: q = (Mat_ScaLAPACK*)Q->data;
99: PetscCallELPARET(elpa_init,20250131);
100: PetscCallELPANOARG(handle = elpa_allocate);
102: /* set parameters of the matrix and its MPI distribution */
103: PetscCallELPA(elpa_set,handle,"na",a->N); /* matrix size */
104: PetscCallELPA(elpa_set,handle,"nev",a->N); /* number of eigenvectors to be computed (1<=nev<=na) */
105: PetscCallELPA(elpa_set,handle,"local_nrows",a->locr); /* number of local rows of the distributed matrix on this MPI task */
106: PetscCallELPA(elpa_set,handle,"local_ncols",a->locc); /* number of local columns of the distributed matrix on this MPI task */
107: PetscCallELPA(elpa_set,handle,"nblk",a->mb); /* size of the BLACS block cyclic distribution */
108: PetscCallELPA(elpa_set,handle,"mpi_comm_parent",MPI_Comm_c2f(PetscObjectComm((PetscObject)eps)));
109: PetscCallELPA(elpa_set,handle,"process_row",a->grid->myrow); /* row coordinate of MPI process */
110: PetscCallELPA(elpa_set,handle,"process_col",a->grid->mycol); /* column coordinate of MPI process */
111: if (B) PetscCallELPA(elpa_set,handle,"blacs_context",a->grid->ictxt);
113: /* setup and set tunable run-time options */
114: PetscCallELPA(elpa_set,handle,"solver",ELPA_SOLVER_2STAGE);
115: /* PetscCallELPA(elpa_print_settings,handle); */
116: PetscCallELPARET(elpa_setup,handle);
118: /* solve the eigenvalue problem */
119: if (B) {
120: b = (Mat_ScaLAPACK*)B->data;
121: PetscCallELPA(elpa_generalized_eigenvectors,handle,a->loc,b->loc,w,q->loc,0);
122: } else PetscCallELPA(elpa_eigenvectors,handle,a->loc,w,q->loc);
124: /* cleanup */
125: PetscCallELPA(elpa_deallocate,handle);
126: PetscCallELPANOARG(elpa_uninit);
128: for (i=0;i<eps->ncv;i++) {
129: eps->eigr[i] = eps->errest[i];
130: eps->errest[i] = PETSC_MACHINE_EPSILON;
131: }
133: PetscCall(BVGetMat(eps->V,&V));
134: PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&V));
135: PetscCall(BVRestoreMat(eps->V,&V));
136: PetscCall(MatDestroy(&Q));
138: eps->nconv = eps->ncv;
139: eps->its = 1;
140: eps->reason = EPS_CONVERGED_TOL;
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: static PetscErrorCode EPSDestroy_ELPA(EPS eps)
145: {
146: PetscFunctionBegin;
147: PetscCall(PetscFree(eps->data));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: static PetscErrorCode EPSReset_ELPA(EPS eps)
152: {
153: EPS_ELPA *ctx = (EPS_ELPA*)eps->data;
155: PetscFunctionBegin;
156: PetscCall(MatDestroy(&ctx->As));
157: PetscCall(MatDestroy(&ctx->Bs));
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: SLEPC_EXTERN PetscErrorCode EPSCreate_ELPA(EPS eps)
162: {
163: EPS_ELPA *ctx;
165: PetscFunctionBegin;
166: PetscCall(PetscNew(&ctx));
167: eps->data = (void*)ctx;
169: eps->categ = EPS_CATEGORY_OTHER;
171: eps->ops->solve = EPSSolve_ELPA;
172: eps->ops->setup = EPSSetUp_ELPA;
173: eps->ops->setupsort = EPSSetUpSort_Basic;
174: eps->ops->destroy = EPSDestroy_ELPA;
175: eps->ops->reset = EPSReset_ELPA;
176: eps->ops->backtransform = EPSBackTransform_Default;
177: eps->ops->setdefaultst = EPSSetDefaultST_NoFactor;
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }