Actual source code: petsc-interface.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }