Line data Source code
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 */
6 :
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"
13 :
14 : static PetscRandom LOBPCG_RandomContext = NULL;
15 :
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_;
20 :
21 : /* type conversion */
22 : n_ = *n;
23 : lda_ = *lda;
24 : info_ = *info;
25 :
26 : LAPACKpotrf_(uplo,&n_,(PetscScalar*)a,&lda_,&info_);
27 :
28 : *info = info_;
29 : return 0;
30 : }
31 :
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_;
35 :
36 : itype_ = *itype;
37 : n_ = *n;
38 : lda_ = *lda;
39 : ldb_ = *ldb;
40 : lwork_ = *lwork;
41 : info_ = *info;
42 :
43 : LAPACKsygv_(&itype_,jobz,uplo,&n_,(PetscScalar*)a,&lda_,(PetscScalar*)b,&ldb_,(PetscScalar*)w,(PetscScalar*)work,&lwork_,&info_);
44 :
45 : *info = info_;
46 : return 0;
47 : }
48 : #else
49 520 : BlopexInt PETSC_zpotrf_interface (char *uplo,BlopexInt *n,komplex *a,BlopexInt* lda,BlopexInt *info)
50 : {
51 520 : PetscBLASInt n_,lda_,info_;
52 :
53 : /* type conversion */
54 520 : n_ = *n;
55 520 : lda_ = (PetscBLASInt)*lda;
56 :
57 520 : LAPACKpotrf_(uplo,&n_,(PetscScalar*)a,&lda_,&info_);
58 :
59 520 : *info = info_;
60 520 : return 0;
61 : }
62 :
63 268 : 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 268 : PetscBLASInt itype_,n_,lda_,ldb_,lwork_,info_;
66 :
67 268 : itype_ = *itype;
68 268 : n_ = *n;
69 268 : lda_ = *lda;
70 268 : ldb_ = *ldb;
71 268 : lwork_ = *lwork;
72 268 : info_ = *info;
73 :
74 268 : LAPACKsygv_(&itype_,jobz,uplo,&n_,(PetscScalar*)a,&lda_,(PetscScalar*)b,&ldb_,(PetscReal*)w,(PetscScalar*)work,&lwork_,(PetscReal*)rwork,&info_);
75 :
76 268 : *info = info_;
77 268 : return 0;
78 : }
79 : #endif
80 :
81 247 : static void *PETSC_MimicVector(void *vvector)
82 : {
83 247 : Vec temp;
84 :
85 247 : PetscCallAbort(PETSC_COMM_SELF,VecDuplicate((Vec)vvector,&temp));
86 247 : return (void*)temp;
87 : }
88 :
89 290 : static BlopexInt PETSC_DestroyVector(void *vvector)
90 : {
91 290 : Vec v = (Vec)vvector;
92 :
93 290 : PetscCall(VecDestroy(&v));
94 : return 0;
95 : }
96 :
97 30457 : static BlopexInt PETSC_InnerProd(void *x,void *y,void *result)
98 : {
99 :
100 30457 : PetscCall(VecDot((Vec)x,(Vec)y,(PetscScalar*)result));
101 : return 0;
102 : }
103 :
104 9257 : static BlopexInt PETSC_CopyVector(void *x,void *y)
105 : {
106 :
107 9257 : PetscCall(VecCopy((Vec)x,(Vec)y));
108 : return 0;
109 : }
110 :
111 12214 : static BlopexInt PETSC_ClearVector(void *x)
112 : {
113 :
114 12214 : PetscCall(VecSet((Vec)x,0.0));
115 : return 0;
116 : }
117 :
118 0 : static BlopexInt PETSC_SetRandomValues(void* v,BlopexInt seed)
119 : {
120 :
121 : /* note: without previous call to LOBPCG_InitRandomContext LOBPCG_RandomContext will be null,
122 : and VecSetRandom will use internal petsc random context */
123 :
124 0 : PetscCall(VecSetRandom((Vec)v,LOBPCG_RandomContext));
125 : return 0;
126 : }
127 :
128 0 : static BlopexInt PETSC_ScaleVector(double alpha,void *x)
129 : {
130 :
131 0 : PetscCall(VecScale((Vec)x,alpha));
132 : return 0;
133 : }
134 :
135 45713 : static BlopexInt PETSC_Axpy(void *alpha,void *x,void *y)
136 : {
137 :
138 45713 : PetscCall(VecAXPY((Vec)y,*(PetscScalar*)alpha,(Vec)x));
139 : return 0;
140 : }
141 :
142 0 : static BlopexInt PETSC_VectorSize(void *x)
143 : {
144 0 : PetscInt N;
145 0 : (void)VecGetSize((Vec)x,&N);
146 0 : return N;
147 : }
148 :
149 6 : int LOBPCG_InitRandomContext(MPI_Comm comm,PetscRandom rand)
150 : {
151 : /* PetscScalar rnd_bound = 1.0; */
152 :
153 6 : if (rand) {
154 0 : PetscCall(PetscObjectReference((PetscObject)rand));
155 0 : PetscCall(PetscRandomDestroy(&LOBPCG_RandomContext));
156 0 : LOBPCG_RandomContext = rand;
157 6 : } else PetscCall(PetscRandomCreate(comm,&LOBPCG_RandomContext));
158 : return 0;
159 : }
160 :
161 6 : int LOBPCG_SetFromOptionsRandomContext(void)
162 : {
163 6 : PetscCall(PetscRandomSetFromOptions(LOBPCG_RandomContext));
164 :
165 : #if defined(PETSC_USE_COMPLEX)
166 6 : 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 : }
172 :
173 6 : int LOBPCG_DestroyRandomContext(void)
174 : {
175 :
176 6 : PetscCall(PetscRandomDestroy(&LOBPCG_RandomContext));
177 : return 0;
178 : }
179 :
180 7 : int PETSCSetupInterpreter(mv_InterfaceInterpreter *i)
181 : {
182 7 : i->CreateVector = PETSC_MimicVector;
183 7 : i->DestroyVector = PETSC_DestroyVector;
184 7 : i->InnerProd = PETSC_InnerProd;
185 7 : i->CopyVector = PETSC_CopyVector;
186 7 : i->ClearVector = PETSC_ClearVector;
187 7 : i->SetRandomValues = PETSC_SetRandomValues;
188 7 : i->ScaleVector = PETSC_ScaleVector;
189 7 : i->Axpy = PETSC_Axpy;
190 7 : i->VectorSize = PETSC_VectorSize;
191 :
192 : /* Multivector part */
193 :
194 7 : i->CreateMultiVector = mv_TempMultiVectorCreateFromSampleVector;
195 7 : i->CopyCreateMultiVector = mv_TempMultiVectorCreateCopy;
196 7 : i->DestroyMultiVector = mv_TempMultiVectorDestroy;
197 :
198 7 : i->Width = mv_TempMultiVectorWidth;
199 7 : i->Height = mv_TempMultiVectorHeight;
200 7 : i->SetMask = mv_TempMultiVectorSetMask;
201 7 : i->CopyMultiVector = mv_TempMultiVectorCopy;
202 7 : i->ClearMultiVector = mv_TempMultiVectorClear;
203 7 : i->SetRandomVectors = mv_TempMultiVectorSetRandom;
204 7 : i->Eval = mv_TempMultiVectorEval;
205 :
206 : #if defined(PETSC_USE_COMPLEX)
207 7 : i->MultiInnerProd = mv_TempMultiVectorByMultiVector_complex;
208 7 : i->MultiInnerProdDiag = mv_TempMultiVectorByMultiVectorDiag_complex;
209 7 : i->MultiVecMat = mv_TempMultiVectorByMatrix_complex;
210 7 : i->MultiVecMatDiag = mv_TempMultiVectorByDiagonal_complex;
211 7 : i->MultiAxpy = mv_TempMultiVectorAxpy_complex;
212 7 : 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
221 :
222 7 : return 0;
223 : }
|
