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