Actual source code: petsc-interface.c
slepc-3.11.2 2019-07-30
1: /* @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ */
2: /* @@@ BLOPEX (version 1.1) LGPL Version 2.1 or above.See www.gnu.org. */
3: /* @@@ Copyright 2010 BLOPEX team http://code.google.com/p/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 "blopex_interpreter.h"
10: #include "blopex_temp_multivector.h"
11: #include "blopex_fortran_matrix.h"
13: static PetscRandom LOBPCG_RandomContext = NULL;
15: BlopexInt PETSC_dpotrf_interface (char *uplo,BlopexInt *n,double *a,BlopexInt * lda,BlopexInt *info)
16: {
17: PetscBLASInt n_,lda_,info_;
19: /* type conversion */
20: n_ = *n;
21: lda_ = *lda;
22: info_ = *info;
24: LAPACKpotrf_(uplo,&n_,(PetscScalar*)a,&lda_,&info_);
26: *info = info_;
27: return 0;
28: }
30: BlopexInt PETSC_zpotrf_interface (char *uplo,BlopexInt *n,komplex *a,BlopexInt* lda,BlopexInt *info)
31: {
32: PetscBLASInt n_,lda_,info_;
34: /* type conversion */
35: n_ = *n;
36: lda_ = (PetscBLASInt)*lda;
38: LAPACKpotrf_(uplo,&n_,(PetscScalar*)a,&lda_,&info_);
40: *info = info_;
41: return 0;
42: }
44: 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)
45: {
46: #if !defined(PETSC_USE_COMPLEX)
47: PetscBLASInt itype_,n_,lda_,ldb_,lwork_,info_;
49: itype_ = *itype;
50: n_ = *n;
51: lda_ = *lda;
52: ldb_ = *ldb;
53: lwork_ = *lwork;
54: info_ = *info;
56: LAPACKsygv_(&itype_,jobz,uplo,&n_,(PetscScalar*)a,&lda_,(PetscScalar*)b,&ldb_,(PetscScalar*)w,(PetscScalar*)work,&lwork_,&info_);
58: *info = info_;
59: #endif
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: #if defined(PETSC_USE_COMPLEX)
66: PetscBLASInt itype_,n_,lda_,ldb_,lwork_,info_;
68: itype_ = *itype;
69: n_ = *n;
70: lda_ = *lda;
71: ldb_ = *ldb;
72: lwork_ = *lwork;
73: info_ = *info;
75: LAPACKsygv_(&itype_,jobz,uplo,&n_,(PetscScalar*)a,&lda_,(PetscScalar*)b,&ldb_,(PetscReal*)w,(PetscScalar*)work,&lwork_,(PetscReal*)rwork,&info_);
77: *info = info_;
78: #endif
79: return 0;
80: }
82: void *PETSC_MimicVector(void *vvector)
83: {
84: PetscErrorCode ierr;
85: Vec temp;
87: VecDuplicate((Vec)vvector,&temp);CHKERRABORT(PETSC_COMM_SELF,ierr);
88: return (void*)temp;
89: }
91: BlopexInt PETSC_DestroyVector(void *vvector)
92: {
94: Vec v = (Vec)vvector;
96: VecDestroy(&v);
97: return 0;
98: }
100: BlopexInt PETSC_InnerProd(void *x,void *y,void *result)
101: {
104: VecDot((Vec)x,(Vec)y,(PetscScalar*)result);
105: return 0;
106: }
108: BlopexInt PETSC_CopyVector(void *x,void *y)
109: {
110: PetscErrorCode ierr;
112: VecCopy((Vec)x,(Vec)y);
113: return 0;
114: }
116: BlopexInt PETSC_ClearVector(void *x)
117: {
118: PetscErrorCode ierr;
120: VecSet((Vec)x,0.0);
121: return 0;
122: }
124: BlopexInt PETSC_SetRandomValues(void* v,BlopexInt seed)
125: {
128: /* note: without previous call to LOBPCG_InitRandomContext LOBPCG_RandomContext will be null,
129: and VecSetRandom will use internal petsc random context */
131: VecSetRandom((Vec)v,LOBPCG_RandomContext);
132: return 0;
133: }
135: BlopexInt PETSC_ScaleVector(double alpha,void *x)
136: {
139: VecScale((Vec)x,alpha);
140: return 0;
141: }
143: BlopexInt PETSC_Axpy(void *alpha,void *x,void *y)
144: {
147: VecAXPY((Vec)y,*(PetscScalar*)alpha,(Vec)x);
148: return 0;
149: }
151: BlopexInt PETSC_VectorSize(void *x)
152: {
153: PetscInt N;
154: VecGetSize((Vec)x,&N);
155: return N;
156: }
158: int LOBPCG_InitRandomContext(MPI_Comm comm,PetscRandom rand)
159: {
161: /* PetscScalar rnd_bound = 1.0; */
163: if (rand) {
164: PetscObjectReference((PetscObject)rand);
165: PetscRandomDestroy(&LOBPCG_RandomContext);
166: LOBPCG_RandomContext = rand;
167: } else {
168: PetscRandomCreate(comm,&LOBPCG_RandomContext);
169: }
170: return 0;
171: }
173: int LOBPCG_SetFromOptionsRandomContext(void)
174: {
176: PetscRandomSetFromOptions(LOBPCG_RandomContext);
178: #if defined(PETSC_USE_COMPLEX)
179: PetscRandomSetInterval(LOBPCG_RandomContext,(PetscScalar)-1.0-1.0*PETSC_i,(PetscScalar)1.0+1.0*PETSC_i);
180: #else
181: PetscRandomSetInterval(LOBPCG_RandomContext,(PetscScalar)-1.0,(PetscScalar)1.0);
182: #endif
183: return 0;
184: }
186: int LOBPCG_DestroyRandomContext(void)
187: {
190: PetscRandomDestroy(&LOBPCG_RandomContext);
191: return 0;
192: }
194: int PETSCSetupInterpreter(mv_InterfaceInterpreter *i)
195: {
196: i->CreateVector = PETSC_MimicVector;
197: i->DestroyVector = PETSC_DestroyVector;
198: i->InnerProd = PETSC_InnerProd;
199: i->CopyVector = PETSC_CopyVector;
200: i->ClearVector = PETSC_ClearVector;
201: i->SetRandomValues = PETSC_SetRandomValues;
202: i->ScaleVector = PETSC_ScaleVector;
203: i->Axpy = PETSC_Axpy;
204: i->VectorSize = PETSC_VectorSize;
206: /* Multivector part */
208: i->CreateMultiVector = mv_TempMultiVectorCreateFromSampleVector;
209: i->CopyCreateMultiVector = mv_TempMultiVectorCreateCopy;
210: i->DestroyMultiVector = mv_TempMultiVectorDestroy;
212: i->Width = mv_TempMultiVectorWidth;
213: i->Height = mv_TempMultiVectorHeight;
214: i->SetMask = mv_TempMultiVectorSetMask;
215: i->CopyMultiVector = mv_TempMultiVectorCopy;
216: i->ClearMultiVector = mv_TempMultiVectorClear;
217: i->SetRandomVectors = mv_TempMultiVectorSetRandom;
218: i->Eval = mv_TempMultiVectorEval;
220: #if defined(PETSC_USE_COMPLEX)
221: i->MultiInnerProd = mv_TempMultiVectorByMultiVector_complex;
222: i->MultiInnerProdDiag = mv_TempMultiVectorByMultiVectorDiag_complex;
223: i->MultiVecMat = mv_TempMultiVectorByMatrix_complex;
224: i->MultiVecMatDiag = mv_TempMultiVectorByDiagonal_complex;
225: i->MultiAxpy = mv_TempMultiVectorAxpy_complex;
226: i->MultiXapy = mv_TempMultiVectorXapy_complex;
227: #else
228: i->MultiInnerProd = mv_TempMultiVectorByMultiVector;
229: i->MultiInnerProdDiag = mv_TempMultiVectorByMultiVectorDiag;
230: i->MultiVecMat = mv_TempMultiVectorByMatrix;
231: i->MultiVecMatDiag = mv_TempMultiVectorByDiagonal;
232: i->MultiAxpy = mv_TempMultiVectorAxpy;
233: i->MultiXapy = mv_TempMultiVectorXapy;
234: #endif
236: return 0;
237: }