Actual source code: elpa.c
slepc-3.21.2 2024-09-25
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: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift));
56: PetscCheck(isshift,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support spectral transformations");
57: eps->ncv = eps->n;
58: if (eps->mpd!=PETSC_DEFAULT) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
59: if (eps->max_it==PETSC_DEFAULT) eps->max_it = 1;
60: if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
61: PetscCheck(eps->which!=EPS_ALL || eps->inta==eps->intb,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support interval computation");
62: EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_STOPPING);
63: EPSCheckIgnored(eps,EPS_FEATURE_EXTRACTION | EPS_FEATURE_CONVERGENCE);
64: PetscCall(EPSAllocateSolution(eps,0));
66: /* convert matrices */
67: PetscCall(MatDestroy(&ctx->As));
68: PetscCall(MatDestroy(&ctx->Bs));
69: PetscCall(STGetNumMatrices(eps->st,&nmat));
70: PetscCall(STGetMatrix(eps->st,0,&A));
71: PetscCall(MatConvert(A,MATSCALAPACK,MAT_INITIAL_MATRIX,&ctx->As));
72: if (nmat>1) {
73: PetscCall(STGetMatrix(eps->st,1,&B));
74: PetscCall(MatConvert(B,MATSCALAPACK,MAT_INITIAL_MATRIX,&ctx->Bs));
75: }
76: PetscCall(STGetShift(eps->st,&shift));
77: if (shift != 0.0) {
78: if (nmat>1) PetscCall(MatAXPY(ctx->As,-shift,ctx->Bs,SAME_NONZERO_PATTERN));
79: else PetscCall(MatShift(ctx->As,-shift));
80: }
81: PetscFunctionReturn(PETSC_SUCCESS);
82: }
84: static PetscErrorCode EPSSolve_ELPA(EPS eps)
85: {
86: EPS_ELPA *ctx = (EPS_ELPA*)eps->data;
87: Mat A = ctx->As,B = ctx->Bs,Q,V;
88: Mat_ScaLAPACK *a = (Mat_ScaLAPACK*)A->data,*b,*q;
89: PetscReal *w = eps->errest; /* used to store real eigenvalues */
90: PetscInt i;
91: elpa_t handle;
93: PetscFunctionBegin;
94: PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Q));
95: q = (Mat_ScaLAPACK*)Q->data;
97: PetscCallELPARET(elpa_init,20200417); /* 20171201 */
98: PetscCallELPANOARG(handle = elpa_allocate);
100: /* set parameters of the matrix and its MPI distribution */
101: PetscCallELPA(elpa_set,handle,"na",a->N); /* matrix size */
102: PetscCallELPA(elpa_set,handle,"nev",a->N); /* number of eigenvectors to be computed (1<=nev<=na) */
103: PetscCallELPA(elpa_set,handle,"local_nrows",a->locr); /* number of local rows of the distributed matrix on this MPI task */
104: PetscCallELPA(elpa_set,handle,"local_ncols",a->locc); /* number of local columns of the distributed matrix on this MPI task */
105: PetscCallELPA(elpa_set,handle,"nblk",a->mb); /* size of the BLACS block cyclic distribution */
106: PetscCallELPA(elpa_set,handle,"mpi_comm_parent",MPI_Comm_c2f(PetscObjectComm((PetscObject)eps)));
107: PetscCallELPA(elpa_set,handle,"process_row",a->grid->myrow); /* row coordinate of MPI process */
108: PetscCallELPA(elpa_set,handle,"process_col",a->grid->mycol); /* column coordinate of MPI process */
109: if (B) PetscCallELPA(elpa_set,handle,"blacs_context",a->grid->ictxt);
111: /* setup and set tunable run-time options */
112: PetscCallELPARET(elpa_setup,handle);
113: PetscCallELPA(elpa_set,handle,"solver",ELPA_SOLVER_2STAGE);
114: /* PetscCallELPA(elpa_print_settings,handle); */
116: /* solve the eigenvalue problem */
117: if (B) {
118: b = (Mat_ScaLAPACK*)B->data;
119: PetscCallELPA(elpa_generalized_eigenvectors,handle,a->loc,b->loc,w,q->loc,0);
120: } else PetscCallELPA(elpa_eigenvectors,handle,a->loc,w,q->loc);
122: /* cleanup */
123: PetscCallELPA(elpa_deallocate,handle);
124: PetscCallELPANOARG(elpa_uninit);
126: for (i=0;i<eps->ncv;i++) {
127: eps->eigr[i] = eps->errest[i];
128: eps->errest[i] = PETSC_MACHINE_EPSILON;
129: }
131: PetscCall(BVGetMat(eps->V,&V));
132: PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&V));
133: PetscCall(BVRestoreMat(eps->V,&V));
134: PetscCall(MatDestroy(&Q));
136: eps->nconv = eps->ncv;
137: eps->its = 1;
138: eps->reason = EPS_CONVERGED_TOL;
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: static PetscErrorCode EPSDestroy_ELPA(EPS eps)
143: {
144: PetscFunctionBegin;
145: PetscCall(PetscFree(eps->data));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: static PetscErrorCode EPSReset_ELPA(EPS eps)
150: {
151: EPS_ELPA *ctx = (EPS_ELPA*)eps->data;
153: PetscFunctionBegin;
154: PetscCall(MatDestroy(&ctx->As));
155: PetscCall(MatDestroy(&ctx->Bs));
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: SLEPC_EXTERN PetscErrorCode EPSCreate_ELPA(EPS eps)
160: {
161: EPS_ELPA *ctx;
163: PetscFunctionBegin;
164: PetscCall(PetscNew(&ctx));
165: eps->data = (void*)ctx;
167: eps->categ = EPS_CATEGORY_OTHER;
169: eps->ops->solve = EPSSolve_ELPA;
170: eps->ops->setup = EPSSetUp_ELPA;
171: eps->ops->setupsort = EPSSetUpSort_Basic;
172: eps->ops->destroy = EPSDestroy_ELPA;
173: eps->ops->reset = EPSReset_ELPA;
174: eps->ops->backtransform = EPSBackTransform_Default;
175: eps->ops->setdefaultst = EPSSetDefaultST_NoFactor;
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }