Actual source code: petsc-interface.c
slepc-3.21.1 2024-04-26
1: /* @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ */
2: /* @@@ BLOPEX (version 1.1) LGPL Version 2.1 or above.See www.gnu.org. */
3: /* @@@ Copyright 2010 BLOPEX team https://github.com/lobpcg/blopex */
4: /* @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ */
5: /* This code was developed by Merico Argentati, Andrew Knyazev, Ilya Lashuk and Evgueni Ovtchinnikov */
7: #include <slepcsys.h>
8: #include <petscblaslapack.h>
9: #include <interpreter.h>
10: #include <temp_multivector.h>
11: #include <fortran_matrix.h>
12: #include "petsc-interface.h"
14: static PetscRandom LOBPCG_RandomContext = NULL;
16: #if !defined(PETSC_USE_COMPLEX)
17: BlopexInt PETSC_dpotrf_interface (char *uplo,BlopexInt *n,double *a,BlopexInt * lda,BlopexInt *info)
18: {
19: PetscBLASInt n_,lda_,info_;
21: /* type conversion */
22: n_ = *n;
23: lda_ = *lda;
24: info_ = *info;
26: LAPACKpotrf_(uplo,&n_,(PetscScalar*)a,&lda_,&info_);
28: *info = info_;
29: return 0;
30: }
32: BlopexInt PETSC_dsygv_interface (BlopexInt *itype,char *jobz,char *uplo,BlopexInt *n,double *a,BlopexInt *lda,double *b,BlopexInt *ldb,double *w,double *work,BlopexInt *lwork,BlopexInt *info)
33: {
34: PetscBLASInt itype_,n_,lda_,ldb_,lwork_,info_;
36: itype_ = *itype;
37: n_ = *n;
38: lda_ = *lda;
39: ldb_ = *ldb;
40: lwork_ = *lwork;
41: info_ = *info;
43: LAPACKsygv_(&itype_,jobz,uplo,&n_,(PetscScalar*)a,&lda_,(PetscScalar*)b,&ldb_,(PetscScalar*)w,(PetscScalar*)work,&lwork_,&info_);
45: *info = info_;
46: return 0;
47: }
48: #else
49: BlopexInt PETSC_zpotrf_interface (char *uplo,BlopexInt *n,komplex *a,BlopexInt* lda,BlopexInt *info)
50: {
51: PetscBLASInt n_,lda_,info_;
53: /* type conversion */
54: n_ = *n;
55: lda_ = (PetscBLASInt)*lda;
57: LAPACKpotrf_(uplo,&n_,(PetscScalar*)a,&lda_,&info_);
59: *info = info_;
60: return 0;
61: }
63: BlopexInt PETSC_zsygv_interface (BlopexInt *itype,char *jobz,char *uplo,BlopexInt *n,komplex *a,BlopexInt *lda,komplex *b,BlopexInt *ldb,double *w,komplex *work,BlopexInt *lwork,double *rwork,BlopexInt *info)
64: {
65: PetscBLASInt itype_,n_,lda_,ldb_,lwork_,info_;
67: itype_ = *itype;
68: n_ = *n;
69: lda_ = *lda;
70: ldb_ = *ldb;
71: lwork_ = *lwork;
72: info_ = *info;
74: LAPACKsygv_(&itype_,jobz,uplo,&n_,(PetscScalar*)a,&lda_,(PetscScalar*)b,&ldb_,(PetscReal*)w,(PetscScalar*)work,&lwork_,(PetscReal*)rwork,&info_);
76: *info = info_;
77: return 0;
78: }
79: #endif
81: static void *PETSC_MimicVector(void *vvector)
82: {
83: Vec temp;
85: PetscCallAbort(PETSC_COMM_SELF,VecDuplicate((Vec)vvector,&temp));
86: return (void*)temp;
87: }
89: static BlopexInt PETSC_DestroyVector(void *vvector)
90: {
91: Vec v = (Vec)vvector;
93: PetscCall(VecDestroy(&v));
94: return 0;
95: }
97: static BlopexInt PETSC_InnerProd(void *x,void *y,void *result)
98: {
100: PetscCall(VecDot((Vec)x,(Vec)y,(PetscScalar*)result));
101: return 0;
102: }
104: static BlopexInt PETSC_CopyVector(void *x,void *y)
105: {
107: PetscCall(VecCopy((Vec)x,(Vec)y));
108: return 0;
109: }
111: static BlopexInt PETSC_ClearVector(void *x)
112: {
114: PetscCall(VecSet((Vec)x,0.0));
115: return 0;
116: }
118: static BlopexInt PETSC_SetRandomValues(void* v,BlopexInt seed)
119: {
121: /* note: without previous call to LOBPCG_InitRandomContext LOBPCG_RandomContext will be null,
122: and VecSetRandom will use internal petsc random context */
124: PetscCall(VecSetRandom((Vec)v,LOBPCG_RandomContext));
125: return 0;
126: }
128: static BlopexInt PETSC_ScaleVector(double alpha,void *x)
129: {
131: PetscCall(VecScale((Vec)x,alpha));
132: return 0;
133: }
135: static BlopexInt PETSC_Axpy(void *alpha,void *x,void *y)
136: {
138: PetscCall(VecAXPY((Vec)y,*(PetscScalar*)alpha,(Vec)x));
139: return 0;
140: }
142: static BlopexInt PETSC_VectorSize(void *x)
143: {
144: PetscInt N;
145: (void)VecGetSize((Vec)x,&N);
146: return N;
147: }
149: int LOBPCG_InitRandomContext(MPI_Comm comm,PetscRandom rand)
150: {
151: /* PetscScalar rnd_bound = 1.0; */
153: if (rand) {
154: PetscCall(PetscObjectReference((PetscObject)rand));
155: PetscCall(PetscRandomDestroy(&LOBPCG_RandomContext));
156: LOBPCG_RandomContext = rand;
157: } else PetscCall(PetscRandomCreate(comm,&LOBPCG_RandomContext));
158: return 0;
159: }
161: int LOBPCG_SetFromOptionsRandomContext(void)
162: {
163: PetscCall(PetscRandomSetFromOptions(LOBPCG_RandomContext));
165: #if defined(PETSC_USE_COMPLEX)
166: PetscCall(PetscRandomSetInterval(LOBPCG_RandomContext,(PetscScalar)PetscCMPLX(-1.0,-1.0),(PetscScalar)PetscCMPLX(1.0,1.0)));
167: #else
168: PetscCall(PetscRandomSetInterval(LOBPCG_RandomContext,(PetscScalar)-1.0,(PetscScalar)1.0));
169: #endif
170: return 0;
171: }
173: int LOBPCG_DestroyRandomContext(void)
174: {
176: PetscCall(PetscRandomDestroy(&LOBPCG_RandomContext));
177: return 0;
178: }
180: int PETSCSetupInterpreter(mv_InterfaceInterpreter *i)
181: {
182: i->CreateVector = PETSC_MimicVector;
183: i->DestroyVector = PETSC_DestroyVector;
184: i->InnerProd = PETSC_InnerProd;
185: i->CopyVector = PETSC_CopyVector;
186: i->ClearVector = PETSC_ClearVector;
187: i->SetRandomValues = PETSC_SetRandomValues;
188: i->ScaleVector = PETSC_ScaleVector;
189: i->Axpy = PETSC_Axpy;
190: i->VectorSize = PETSC_VectorSize;
192: /* Multivector part */
194: i->CreateMultiVector = mv_TempMultiVectorCreateFromSampleVector;
195: i->CopyCreateMultiVector = mv_TempMultiVectorCreateCopy;
196: i->DestroyMultiVector = mv_TempMultiVectorDestroy;
198: i->Width = mv_TempMultiVectorWidth;
199: i->Height = mv_TempMultiVectorHeight;
200: i->SetMask = mv_TempMultiVectorSetMask;
201: i->CopyMultiVector = mv_TempMultiVectorCopy;
202: i->ClearMultiVector = mv_TempMultiVectorClear;
203: i->SetRandomVectors = mv_TempMultiVectorSetRandom;
204: i->Eval = mv_TempMultiVectorEval;
206: #if defined(PETSC_USE_COMPLEX)
207: i->MultiInnerProd = mv_TempMultiVectorByMultiVector_complex;
208: i->MultiInnerProdDiag = mv_TempMultiVectorByMultiVectorDiag_complex;
209: i->MultiVecMat = mv_TempMultiVectorByMatrix_complex;
210: i->MultiVecMatDiag = mv_TempMultiVectorByDiagonal_complex;
211: i->MultiAxpy = mv_TempMultiVectorAxpy_complex;
212: i->MultiXapy = mv_TempMultiVectorXapy_complex;
213: #else
214: i->MultiInnerProd = mv_TempMultiVectorByMultiVector;
215: i->MultiInnerProdDiag = mv_TempMultiVectorByMultiVectorDiag;
216: i->MultiVecMat = mv_TempMultiVectorByMatrix;
217: i->MultiVecMatDiag = mv_TempMultiVectorByDiagonal;
218: i->MultiAxpy = mv_TempMultiVectorAxpy;
219: i->MultiXapy = mv_TempMultiVectorXapy;
220: #endif
222: return 0;
223: }