Actual source code: petsc-interface.c
slepc-3.18.0 2022-10-01
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 <petscvec.h>
8: #include <petscblaslapack.h>
9: #include <interpreter.h>
10: #include <temp_multivector.h>
11: #include <fortran_matrix.h>
13: static PetscRandom LOBPCG_RandomContext = NULL;
15: #if !defined(PETSC_USE_COMPLEX)
16: BlopexInt PETSC_dpotrf_interface (char *uplo,BlopexInt *n,double *a,BlopexInt * lda,BlopexInt *info)
17: {
18: PetscBLASInt n_,lda_,info_;
20: /* type conversion */
21: n_ = *n;
22: lda_ = *lda;
23: info_ = *info;
25: LAPACKpotrf_(uplo,&n_,(PetscScalar*)a,&lda_,&info_);
27: *info = info_;
28: return 0;
29: }
31: 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)
32: {
33: PetscBLASInt itype_,n_,lda_,ldb_,lwork_,info_;
35: itype_ = *itype;
36: n_ = *n;
37: lda_ = *lda;
38: ldb_ = *ldb;
39: lwork_ = *lwork;
40: info_ = *info;
42: LAPACKsygv_(&itype_,jobz,uplo,&n_,(PetscScalar*)a,&lda_,(PetscScalar*)b,&ldb_,(PetscScalar*)w,(PetscScalar*)work,&lwork_,&info_);
44: *info = info_;
45: return 0;
46: }
47: #else
48: BlopexInt PETSC_zpotrf_interface (char *uplo,BlopexInt *n,komplex *a,BlopexInt* lda,BlopexInt *info)
49: {
50: PetscBLASInt n_,lda_,info_;
52: /* type conversion */
53: n_ = *n;
54: lda_ = (PetscBLASInt)*lda;
56: LAPACKpotrf_(uplo,&n_,(PetscScalar*)a,&lda_,&info_);
58: *info = info_;
59: return 0;
60: }
62: 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)
63: {
64: PetscBLASInt itype_,n_,lda_,ldb_,lwork_,info_;
66: itype_ = *itype;
67: n_ = *n;
68: lda_ = *lda;
69: ldb_ = *ldb;
70: lwork_ = *lwork;
71: info_ = *info;
73: LAPACKsygv_(&itype_,jobz,uplo,&n_,(PetscScalar*)a,&lda_,(PetscScalar*)b,&ldb_,(PetscReal*)w,(PetscScalar*)work,&lwork_,(PetscReal*)rwork,&info_);
75: *info = info_;
76: return 0;
77: }
78: #endif
80: void *PETSC_MimicVector(void *vvector)
81: {
82: Vec temp;
84: PETSC_COMM_SELF,VecDuplicate((Vec)vvector,&temp);
85: return (void*)temp;
86: }
88: BlopexInt PETSC_DestroyVector(void *vvector)
89: {
90: Vec v = (Vec)vvector;
92: VecDestroy(&v);
93: return 0;
94: }
96: BlopexInt PETSC_InnerProd(void *x,void *y,void *result)
97: {
99: VecDot((Vec)x,(Vec)y,(PetscScalar*)result);
100: return 0;
101: }
103: BlopexInt PETSC_CopyVector(void *x,void *y)
104: {
106: VecCopy((Vec)x,(Vec)y);
107: return 0;
108: }
110: BlopexInt PETSC_ClearVector(void *x)
111: {
113: VecSet((Vec)x,0.0);
114: return 0;
115: }
117: BlopexInt PETSC_SetRandomValues(void* v,BlopexInt seed)
118: {
120: /* note: without previous call to LOBPCG_InitRandomContext LOBPCG_RandomContext will be null,
121: and VecSetRandom will use internal petsc random context */
123: VecSetRandom((Vec)v,LOBPCG_RandomContext);
124: return 0;
125: }
127: BlopexInt PETSC_ScaleVector(double alpha,void *x)
128: {
130: VecScale((Vec)x,alpha);
131: return 0;
132: }
134: BlopexInt PETSC_Axpy(void *alpha,void *x,void *y)
135: {
137: VecAXPY((Vec)y,*(PetscScalar*)alpha,(Vec)x);
138: return 0;
139: }
141: BlopexInt PETSC_VectorSize(void *x)
142: {
143: PetscInt N;
144: VecGetSize((Vec)x,&N);
145: return N;
146: }
148: int LOBPCG_InitRandomContext(MPI_Comm comm,PetscRandom rand)
149: {
150: /* PetscScalar rnd_bound = 1.0; */
152: if (rand) {
153: PetscObjectReference((PetscObject)rand);
154: PetscRandomDestroy(&LOBPCG_RandomContext);
155: LOBPCG_RandomContext = rand;
156: } else PetscRandomCreate(comm,&LOBPCG_RandomContext);
157: return 0;
158: }
160: int LOBPCG_SetFromOptionsRandomContext(void)
161: {
162: PetscRandomSetFromOptions(LOBPCG_RandomContext);
164: #if defined(PETSC_USE_COMPLEX)
165: PetscRandomSetInterval(LOBPCG_RandomContext,(PetscScalar)PetscCMPLX(-1.0,-1.0),(PetscScalar)PetscCMPLX(1.0,1.0));
166: #else
167: PetscRandomSetInterval(LOBPCG_RandomContext,(PetscScalar)-1.0,(PetscScalar)1.0);
168: #endif
169: return 0;
170: }
172: int LOBPCG_DestroyRandomContext(void)
173: {
175: PetscRandomDestroy(&LOBPCG_RandomContext);
176: return 0;
177: }
179: int PETSCSetupInterpreter(mv_InterfaceInterpreter *i)
180: {
181: i->CreateVector = PETSC_MimicVector;
182: i->DestroyVector = PETSC_DestroyVector;
183: i->InnerProd = PETSC_InnerProd;
184: i->CopyVector = PETSC_CopyVector;
185: i->ClearVector = PETSC_ClearVector;
186: i->SetRandomValues = PETSC_SetRandomValues;
187: i->ScaleVector = PETSC_ScaleVector;
188: i->Axpy = PETSC_Axpy;
189: i->VectorSize = PETSC_VectorSize;
191: /* Multivector part */
193: i->CreateMultiVector = mv_TempMultiVectorCreateFromSampleVector;
194: i->CopyCreateMultiVector = mv_TempMultiVectorCreateCopy;
195: i->DestroyMultiVector = mv_TempMultiVectorDestroy;
197: i->Width = mv_TempMultiVectorWidth;
198: i->Height = mv_TempMultiVectorHeight;
199: i->SetMask = mv_TempMultiVectorSetMask;
200: i->CopyMultiVector = mv_TempMultiVectorCopy;
201: i->ClearMultiVector = mv_TempMultiVectorClear;
202: i->SetRandomVectors = mv_TempMultiVectorSetRandom;
203: i->Eval = mv_TempMultiVectorEval;
205: #if defined(PETSC_USE_COMPLEX)
206: i->MultiInnerProd = mv_TempMultiVectorByMultiVector_complex;
207: i->MultiInnerProdDiag = mv_TempMultiVectorByMultiVectorDiag_complex;
208: i->MultiVecMat = mv_TempMultiVectorByMatrix_complex;
209: i->MultiVecMatDiag = mv_TempMultiVectorByDiagonal_complex;
210: i->MultiAxpy = mv_TempMultiVectorAxpy_complex;
211: i->MultiXapy = mv_TempMultiVectorXapy_complex;
212: #else
213: i->MultiInnerProd = mv_TempMultiVectorByMultiVector;
214: i->MultiInnerProdDiag = mv_TempMultiVectorByMultiVectorDiag;
215: i->MultiVecMat = mv_TempMultiVectorByMatrix;
216: i->MultiVecMatDiag = mv_TempMultiVectorByDiagonal;
217: i->MultiAxpy = mv_TempMultiVectorAxpy;
218: i->MultiXapy = mv_TempMultiVectorXapy;
219: #endif
221: return 0;
222: }