Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 : /*
11 : BV implemented as a single Vec
12 : */
13 :
14 : #include <slepc/private/bvimpl.h>
15 : #include "svec.h"
16 :
17 100 : static PetscErrorCode BVMult_Svec(BV Y,PetscScalar alpha,PetscScalar beta,BV X,Mat Q)
18 : {
19 100 : BV_SVEC *y = (BV_SVEC*)Y->data,*x = (BV_SVEC*)X->data;
20 100 : const PetscScalar *px,*q;
21 100 : PetscScalar *py;
22 100 : PetscInt ldq;
23 :
24 100 : PetscFunctionBegin;
25 100 : PetscCall(VecGetArrayRead(x->v,&px));
26 100 : PetscCall(VecGetArray(y->v,&py));
27 100 : if (Q) {
28 80 : PetscCall(MatDenseGetLDA(Q,&ldq));
29 80 : PetscCall(MatDenseGetArrayRead(Q,&q));
30 80 : PetscCall(BVMult_BLAS_Private(Y,Y->n,Y->k-Y->l,X->k-X->l,alpha,px+(X->nc+X->l)*X->ld,X->ld,q+Y->l*ldq+X->l,ldq,beta,py+(Y->nc+Y->l)*Y->ld,Y->ld));
31 80 : PetscCall(MatDenseRestoreArrayRead(Q,&q));
32 20 : } else PetscCall(BVAXPY_BLAS_Private(Y,Y->n,Y->k-Y->l,alpha,px+(X->nc+X->l)*X->ld,X->ld,beta,py+(Y->nc+Y->l)*Y->ld,Y->ld));
33 100 : PetscCall(VecRestoreArrayRead(x->v,&px));
34 100 : PetscCall(VecRestoreArray(y->v,&py));
35 100 : PetscFunctionReturn(PETSC_SUCCESS);
36 : }
37 :
38 347 : static PetscErrorCode BVMultVec_Svec(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar *q)
39 : {
40 347 : BV_SVEC *x = (BV_SVEC*)X->data;
41 347 : PetscScalar *px,*py,*qq=q;
42 :
43 347 : PetscFunctionBegin;
44 347 : PetscCall(VecGetArray(x->v,&px));
45 347 : PetscCall(VecGetArray(y,&py));
46 347 : if (!q) PetscCall(VecGetArray(X->buffer,&qq));
47 347 : PetscCall(BVMultVec_BLAS_Private(X,X->n,X->k-X->l,alpha,px+(X->nc+X->l)*X->ld,X->ld,qq,beta,py));
48 347 : if (!q) PetscCall(VecRestoreArray(X->buffer,&qq));
49 347 : PetscCall(VecRestoreArray(x->v,&px));
50 347 : PetscCall(VecRestoreArray(y,&py));
51 347 : PetscFunctionReturn(PETSC_SUCCESS);
52 : }
53 :
54 46 : static PetscErrorCode BVMultInPlace_Svec(BV V,Mat Q,PetscInt s,PetscInt e)
55 : {
56 46 : BV_SVEC *ctx = (BV_SVEC*)V->data;
57 46 : PetscScalar *pv;
58 46 : const PetscScalar *q;
59 46 : PetscInt ldq;
60 :
61 46 : PetscFunctionBegin;
62 46 : if (s>=e || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
63 46 : PetscCall(MatDenseGetLDA(Q,&ldq));
64 46 : PetscCall(VecGetArray(ctx->v,&pv));
65 46 : PetscCall(MatDenseGetArrayRead(Q,&q));
66 46 : PetscCall(BVMultInPlace_BLAS_Private(V,V->n,V->k-V->l,s-V->l,e-V->l,pv+(V->nc+V->l)*V->ld,V->ld,q+V->l*ldq+V->l,ldq,PETSC_FALSE));
67 46 : PetscCall(MatDenseRestoreArrayRead(Q,&q));
68 46 : PetscCall(VecRestoreArray(ctx->v,&pv));
69 46 : PetscFunctionReturn(PETSC_SUCCESS);
70 : }
71 :
72 1 : static PetscErrorCode BVMultInPlaceHermitianTranspose_Svec(BV V,Mat Q,PetscInt s,PetscInt e)
73 : {
74 1 : BV_SVEC *ctx = (BV_SVEC*)V->data;
75 1 : PetscScalar *pv;
76 1 : const PetscScalar *q;
77 1 : PetscInt ldq;
78 :
79 1 : PetscFunctionBegin;
80 1 : if (s>=e || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
81 1 : PetscCall(MatDenseGetLDA(Q,&ldq));
82 1 : PetscCall(VecGetArray(ctx->v,&pv));
83 1 : PetscCall(MatDenseGetArrayRead(Q,&q));
84 1 : PetscCall(BVMultInPlace_BLAS_Private(V,V->n,V->k-V->l,s-V->l,e-V->l,pv+(V->nc+V->l)*V->ld,V->ld,q+V->l*ldq+V->l,ldq,PETSC_TRUE));
85 1 : PetscCall(MatDenseRestoreArrayRead(Q,&q));
86 1 : PetscCall(VecRestoreArray(ctx->v,&pv));
87 1 : PetscFunctionReturn(PETSC_SUCCESS);
88 : }
89 :
90 243 : static PetscErrorCode BVDot_Svec(BV X,BV Y,Mat M)
91 : {
92 243 : BV_SVEC *x = (BV_SVEC*)X->data,*y = (BV_SVEC*)Y->data;
93 243 : const PetscScalar *px,*py;
94 243 : PetscScalar *m;
95 243 : PetscInt ldm;
96 :
97 243 : PetscFunctionBegin;
98 243 : PetscCall(MatDenseGetLDA(M,&ldm));
99 243 : PetscCall(VecGetArrayRead(x->v,&px));
100 243 : PetscCall(VecGetArrayRead(y->v,&py));
101 243 : PetscCall(MatDenseGetArray(M,&m));
102 243 : PetscCall(BVDot_BLAS_Private(X,Y->k-Y->l,X->k-X->l,X->n,py+(Y->nc+Y->l)*Y->ld,Y->ld,px+(X->nc+X->l)*X->ld,X->ld,m+X->l*ldm+Y->l,ldm,x->mpi));
103 243 : PetscCall(MatDenseRestoreArray(M,&m));
104 243 : PetscCall(VecRestoreArrayRead(x->v,&px));
105 243 : PetscCall(VecRestoreArrayRead(y->v,&py));
106 243 : PetscFunctionReturn(PETSC_SUCCESS);
107 : }
108 :
109 405 : static PetscErrorCode BVDotVec_Svec(BV X,Vec y,PetscScalar *q)
110 : {
111 405 : BV_SVEC *x = (BV_SVEC*)X->data;
112 405 : const PetscScalar *px,*py;
113 405 : PetscScalar *qq=q;
114 405 : Vec z = y;
115 :
116 405 : PetscFunctionBegin;
117 405 : if (PetscUnlikely(X->matrix)) {
118 62 : PetscCall(BV_IPMatMult(X,y));
119 62 : z = X->Bx;
120 : }
121 405 : PetscCall(VecGetArrayRead(x->v,&px));
122 405 : PetscCall(VecGetArrayRead(z,&py));
123 405 : if (!q) PetscCall(VecGetArray(X->buffer,&qq));
124 405 : PetscCall(BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->ld,X->ld,py,qq,x->mpi));
125 405 : if (!q) PetscCall(VecRestoreArray(X->buffer,&qq));
126 405 : PetscCall(VecRestoreArrayRead(z,&py));
127 405 : PetscCall(VecRestoreArrayRead(x->v,&px));
128 405 : PetscFunctionReturn(PETSC_SUCCESS);
129 : }
130 :
131 4 : static PetscErrorCode BVDotVec_Local_Svec(BV X,Vec y,PetscScalar *m)
132 : {
133 4 : BV_SVEC *x = (BV_SVEC*)X->data;
134 4 : PetscScalar *px,*py;
135 4 : Vec z = y;
136 :
137 4 : PetscFunctionBegin;
138 4 : if (PetscUnlikely(X->matrix)) {
139 0 : PetscCall(BV_IPMatMult(X,y));
140 0 : z = X->Bx;
141 : }
142 4 : PetscCall(VecGetArray(x->v,&px));
143 4 : PetscCall(VecGetArray(z,&py));
144 4 : PetscCall(BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->ld,X->ld,py,m,PETSC_FALSE));
145 4 : PetscCall(VecRestoreArray(z,&py));
146 4 : PetscCall(VecRestoreArray(x->v,&px));
147 4 : PetscFunctionReturn(PETSC_SUCCESS);
148 : }
149 :
150 328 : static PetscErrorCode BVScale_Svec(BV bv,PetscInt j,PetscScalar alpha)
151 : {
152 328 : BV_SVEC *ctx = (BV_SVEC*)bv->data;
153 328 : PetscScalar *array;
154 :
155 328 : PetscFunctionBegin;
156 328 : if (!bv->n) PetscFunctionReturn(PETSC_SUCCESS);
157 328 : PetscCall(VecGetArray(ctx->v,&array));
158 328 : if (PetscUnlikely(j<0)) PetscCall(BVScale_BLAS_Private(bv,(bv->k-bv->l)*bv->ld,array+(bv->nc+bv->l)*bv->ld,alpha));
159 324 : else PetscCall(BVScale_BLAS_Private(bv,bv->n,array+(bv->nc+j)*bv->ld,alpha));
160 328 : PetscCall(VecRestoreArray(ctx->v,&array));
161 328 : PetscFunctionReturn(PETSC_SUCCESS);
162 : }
163 :
164 112 : static PetscErrorCode BVNorm_Svec(BV bv,PetscInt j,NormType type,PetscReal *val)
165 : {
166 112 : BV_SVEC *ctx = (BV_SVEC*)bv->data;
167 112 : const PetscScalar *array;
168 :
169 112 : PetscFunctionBegin;
170 112 : PetscCall(VecGetArrayRead(ctx->v,&array));
171 112 : if (PetscUnlikely(j<0)) PetscCall(BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->ld,bv->ld,type,val,ctx->mpi));
172 42 : else PetscCall(BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->ld,bv->ld,type,val,ctx->mpi));
173 112 : PetscCall(VecRestoreArrayRead(ctx->v,&array));
174 112 : PetscFunctionReturn(PETSC_SUCCESS);
175 : }
176 :
177 2 : static PetscErrorCode BVNorm_Local_Svec(BV bv,PetscInt j,NormType type,PetscReal *val)
178 : {
179 2 : BV_SVEC *ctx = (BV_SVEC*)bv->data;
180 2 : const PetscScalar *array;
181 :
182 2 : PetscFunctionBegin;
183 2 : PetscCall(VecGetArrayRead(ctx->v,&array));
184 2 : if (PetscUnlikely(j<0)) PetscCall(BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->ld,bv->ld,type,val,PETSC_FALSE));
185 2 : else PetscCall(BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->ld,bv->ld,type,val,PETSC_FALSE));
186 2 : PetscCall(VecRestoreArrayRead(ctx->v,&array));
187 2 : PetscFunctionReturn(PETSC_SUCCESS);
188 : }
189 :
190 6 : static PetscErrorCode BVNormalize_Svec(BV bv,PetscScalar *eigi)
191 : {
192 6 : BV_SVEC *ctx = (BV_SVEC*)bv->data;
193 6 : PetscScalar *array,*wi=NULL;
194 :
195 6 : PetscFunctionBegin;
196 6 : PetscCall(VecGetArray(ctx->v,&array));
197 6 : if (eigi) wi = eigi+bv->l;
198 6 : PetscCall(BVNormalize_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->ld,bv->ld,wi,ctx->mpi));
199 6 : PetscCall(VecRestoreArray(ctx->v,&array));
200 6 : PetscFunctionReturn(PETSC_SUCCESS);
201 : }
202 :
203 111 : static PetscErrorCode BVMatMult_Svec(BV V,Mat A,BV W)
204 : {
205 111 : PetscInt j;
206 111 : Mat Vmat,Wmat;
207 111 : Vec vv,ww;
208 :
209 111 : PetscFunctionBegin;
210 111 : if (V->vmm) {
211 100 : PetscCall(BVGetMat(V,&Vmat));
212 100 : PetscCall(BVGetMat(W,&Wmat));
213 100 : PetscCall(MatProductCreateWithMat(A,Vmat,NULL,Wmat));
214 100 : PetscCall(MatProductSetType(Wmat,MATPRODUCT_AB));
215 100 : PetscCall(MatProductSetFromOptions(Wmat));
216 100 : PetscCall(MatProductSymbolic(Wmat));
217 100 : PetscCall(MatProductNumeric(Wmat));
218 100 : PetscCall(MatProductClear(Wmat));
219 100 : PetscCall(BVRestoreMat(V,&Vmat));
220 100 : PetscCall(BVRestoreMat(W,&Wmat));
221 : } else {
222 77 : for (j=0;j<V->k-V->l;j++) {
223 66 : PetscCall(BVGetColumn(V,V->l+j,&vv));
224 66 : PetscCall(BVGetColumn(W,W->l+j,&ww));
225 66 : PetscCall(MatMult(A,vv,ww));
226 66 : PetscCall(BVRestoreColumn(V,V->l+j,&vv));
227 66 : PetscCall(BVRestoreColumn(W,W->l+j,&ww));
228 : }
229 : }
230 111 : PetscFunctionReturn(PETSC_SUCCESS);
231 : }
232 :
233 99 : static PetscErrorCode BVCopy_Svec(BV V,BV W)
234 : {
235 99 : BV_SVEC *v = (BV_SVEC*)V->data,*w = (BV_SVEC*)W->data;
236 99 : const PetscScalar *pv;
237 99 : PetscScalar *pw;
238 99 : PetscInt j;
239 :
240 99 : PetscFunctionBegin;
241 99 : PetscCall(VecGetArrayRead(v->v,&pv));
242 99 : PetscCall(VecGetArray(w->v,&pw));
243 742 : for (j=0;j<V->k-V->l;j++) PetscCall(PetscArraycpy(pw+(W->nc+W->l+j)*W->ld,pv+(V->nc+V->l+j)*V->ld,V->n));
244 99 : PetscCall(VecRestoreArrayRead(v->v,&pv));
245 99 : PetscCall(VecRestoreArray(w->v,&pw));
246 99 : PetscFunctionReturn(PETSC_SUCCESS);
247 : }
248 :
249 2 : static PetscErrorCode BVCopyColumn_Svec(BV V,PetscInt j,PetscInt i)
250 : {
251 2 : BV_SVEC *v = (BV_SVEC*)V->data;
252 2 : PetscScalar *pv;
253 :
254 2 : PetscFunctionBegin;
255 2 : PetscCall(VecGetArray(v->v,&pv));
256 2 : PetscCall(PetscArraycpy(pv+(V->nc+i)*V->ld,pv+(V->nc+j)*V->ld,V->n));
257 2 : PetscCall(VecRestoreArray(v->v,&pv));
258 2 : PetscFunctionReturn(PETSC_SUCCESS);
259 : }
260 :
261 8 : static PetscErrorCode BVResize_Svec(BV bv,PetscInt m,PetscBool copy)
262 : {
263 8 : BV_SVEC *ctx = (BV_SVEC*)bv->data;
264 8 : PetscScalar *pnew;
265 8 : const PetscScalar *pv;
266 8 : PetscInt bs;
267 8 : Vec vnew;
268 8 : char str[50];
269 :
270 8 : PetscFunctionBegin;
271 8 : PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
272 8 : PetscCall(VecCreate(PetscObjectComm((PetscObject)bv),&vnew));
273 8 : PetscCall(VecSetType(vnew,bv->vtype));
274 8 : PetscCall(VecSetSizes(vnew,m*bv->ld,PETSC_DECIDE));
275 8 : PetscCall(VecSetBlockSize(vnew,bs));
276 8 : if (((PetscObject)bv)->name) {
277 8 : PetscCall(PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name));
278 8 : PetscCall(PetscObjectSetName((PetscObject)vnew,str));
279 : }
280 8 : if (copy) {
281 2 : PetscCall(VecGetArrayRead(ctx->v,&pv));
282 2 : PetscCall(VecGetArray(vnew,&pnew));
283 2 : PetscCall(PetscArraycpy(pnew,pv,PetscMin(m,bv->m)*bv->ld));
284 2 : PetscCall(VecRestoreArrayRead(ctx->v,&pv));
285 2 : PetscCall(VecRestoreArray(vnew,&pnew));
286 : }
287 8 : PetscCall(VecDestroy(&ctx->v));
288 8 : ctx->v = vnew;
289 8 : PetscFunctionReturn(PETSC_SUCCESS);
290 : }
291 :
292 2761 : static PetscErrorCode BVGetColumn_Svec(BV bv,PetscInt j,Vec *v)
293 : {
294 2761 : BV_SVEC *ctx = (BV_SVEC*)bv->data;
295 2761 : PetscScalar *pv;
296 2761 : PetscInt l;
297 :
298 2761 : PetscFunctionBegin;
299 2761 : l = BVAvailableVec;
300 2761 : PetscCall(VecGetArray(ctx->v,&pv));
301 2761 : PetscCall(VecPlaceArray(bv->cv[l],pv+(bv->nc+j)*bv->ld));
302 2761 : PetscFunctionReturn(PETSC_SUCCESS);
303 : }
304 :
305 2761 : static PetscErrorCode BVRestoreColumn_Svec(BV bv,PetscInt j,Vec *v)
306 : {
307 2761 : BV_SVEC *ctx = (BV_SVEC*)bv->data;
308 2761 : PetscInt l;
309 :
310 2761 : PetscFunctionBegin;
311 2761 : l = (j==bv->ci[0])? 0: 1;
312 2761 : PetscCall(VecResetArray(bv->cv[l]));
313 2761 : PetscCall(VecRestoreArray(ctx->v,NULL));
314 2761 : PetscFunctionReturn(PETSC_SUCCESS);
315 : }
316 :
317 253 : static PetscErrorCode BVGetArray_Svec(BV bv,PetscScalar **a)
318 : {
319 253 : BV_SVEC *ctx = (BV_SVEC*)bv->data;
320 :
321 253 : PetscFunctionBegin;
322 253 : PetscCall(VecGetArray(ctx->v,a));
323 253 : PetscFunctionReturn(PETSC_SUCCESS);
324 : }
325 :
326 253 : static PetscErrorCode BVRestoreArray_Svec(BV bv,PetscScalar **a)
327 : {
328 253 : BV_SVEC *ctx = (BV_SVEC*)bv->data;
329 :
330 253 : PetscFunctionBegin;
331 253 : PetscCall(VecRestoreArray(ctx->v,a));
332 253 : PetscFunctionReturn(PETSC_SUCCESS);
333 : }
334 :
335 8 : static PetscErrorCode BVGetArrayRead_Svec(BV bv,const PetscScalar **a)
336 : {
337 8 : BV_SVEC *ctx = (BV_SVEC*)bv->data;
338 :
339 8 : PetscFunctionBegin;
340 8 : PetscCall(VecGetArrayRead(ctx->v,a));
341 8 : PetscFunctionReturn(PETSC_SUCCESS);
342 : }
343 :
344 8 : static PetscErrorCode BVRestoreArrayRead_Svec(BV bv,const PetscScalar **a)
345 : {
346 8 : BV_SVEC *ctx = (BV_SVEC*)bv->data;
347 :
348 8 : PetscFunctionBegin;
349 8 : PetscCall(VecRestoreArrayRead(ctx->v,a));
350 8 : PetscFunctionReturn(PETSC_SUCCESS);
351 : }
352 :
353 12 : static PetscErrorCode BVView_Svec(BV bv,PetscViewer viewer)
354 : {
355 12 : PetscInt j;
356 12 : Vec v;
357 12 : PetscViewerFormat format;
358 12 : PetscBool isascii,ismatlab=PETSC_FALSE;
359 12 : const char *bvname,*name;
360 :
361 12 : PetscFunctionBegin;
362 12 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
363 12 : if (isascii) {
364 12 : PetscCall(PetscViewerGetFormat(viewer,&format));
365 12 : if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
366 10 : if (format == PETSC_VIEWER_ASCII_MATLAB) ismatlab = PETSC_TRUE;
367 : }
368 0 : if (ismatlab) {
369 0 : PetscCall(PetscObjectGetName((PetscObject)bv,&bvname));
370 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s=[];\n",bvname));
371 : }
372 48 : for (j=0;j<bv->m;j++) {
373 38 : PetscCall(BVGetColumn(bv,j,&v));
374 38 : PetscCall(VecView(v,viewer));
375 38 : if (ismatlab) {
376 0 : PetscCall(PetscObjectGetName((PetscObject)v,&name));
377 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s=[%s,%s];clear %s\n",bvname,bvname,name,name));
378 : }
379 38 : PetscCall(BVRestoreColumn(bv,j,&v));
380 : }
381 10 : PetscFunctionReturn(PETSC_SUCCESS);
382 : }
383 :
384 348 : static PetscErrorCode BVDestroy_Svec(BV bv)
385 : {
386 348 : BV_SVEC *ctx = (BV_SVEC*)bv->data;
387 :
388 348 : PetscFunctionBegin;
389 348 : PetscCall(VecDestroy(&ctx->v));
390 348 : PetscCall(VecDestroy(&bv->cv[0]));
391 348 : PetscCall(VecDestroy(&bv->cv[1]));
392 348 : PetscCall(PetscFree(bv->data));
393 348 : bv->cuda = PETSC_FALSE;
394 348 : PetscFunctionReturn(PETSC_SUCCESS);
395 : }
396 :
397 348 : SLEPC_EXTERN PetscErrorCode BVCreate_Svec(BV bv)
398 : {
399 348 : BV_SVEC *ctx;
400 348 : PetscInt nloc,N,bs,tglobal=0,tlocal,lsplit,j,lda;
401 348 : PetscBool seq,isdense;
402 348 : PetscScalar *vv;
403 348 : const PetscScalar *aa,*array,*ptr;
404 348 : char str[50];
405 348 : BV parent;
406 348 : Vec vpar;
407 348 : MatType mtype;
408 :
409 348 : PetscFunctionBegin;
410 348 : PetscCall(PetscNew(&ctx));
411 348 : bv->data = (void*)ctx;
412 :
413 348 : PetscCall(PetscStrcmpAny(bv->vtype,&bv->cuda,VECSEQCUDA,VECMPICUDA,""));
414 348 : PetscCall(PetscStrcmpAny(bv->vtype,&ctx->mpi,VECMPI,VECMPICUDA,""));
415 :
416 348 : PetscCall(PetscStrcmp(bv->vtype,VECSEQ,&seq));
417 348 : PetscCheck(seq || ctx->mpi || bv->cuda,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"BVSVEC does not support the requested vector type: %s",bv->vtype);
418 :
419 348 : PetscCall(PetscLayoutGetLocalSize(bv->map,&nloc));
420 348 : PetscCall(PetscLayoutGetSize(bv->map,&N));
421 348 : PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
422 348 : PetscCall(BV_SetDefaultLD(bv,nloc));
423 348 : tlocal = bv->m*bv->ld;
424 348 : PetscCall(PetscIntMultError(bv->m,N,&tglobal)); /* just to check integer overflow */
425 :
426 348 : if (PetscUnlikely(bv->issplit)) {
427 : /* split BV: create Vec sharing the memory of the parent BV */
428 80 : parent = bv->splitparent;
429 80 : lsplit = parent->lsplit;
430 80 : vpar = ((BV_SVEC*)parent->data)->v;
431 80 : if (bv->cuda) {
432 : #if defined(PETSC_HAVE_CUDA)
433 : PetscCall(VecCUDAGetArrayRead(vpar,&array));
434 : ptr = (bv->issplit==1)? array: array+lsplit*bv->ld;
435 : PetscCall(VecCUDARestoreArrayRead(vpar,&array));
436 : if (ctx->mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)bv),bs,tlocal,PETSC_DECIDE,NULL,&ctx->v));
437 : else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)bv),bs,tlocal,NULL,&ctx->v));
438 : PetscCall(VecCUDAPlaceArray(ctx->v,ptr));
439 : #endif
440 : } else {
441 80 : PetscCall(VecGetArrayRead(vpar,&array));
442 80 : ptr = (bv->issplit==1)? array: array+lsplit*bv->ld;
443 80 : PetscCall(VecRestoreArrayRead(vpar,&array));
444 80 : if (ctx->mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv),bs,tlocal,PETSC_DECIDE,NULL,&ctx->v));
445 0 : else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv),bs,tlocal,NULL,&ctx->v));
446 80 : PetscCall(VecPlaceArray(ctx->v,ptr));
447 : }
448 : } else {
449 : /* regular BV: create Vec to store the BV entries */
450 268 : PetscCall(VecCreate(PetscObjectComm((PetscObject)bv),&ctx->v));
451 268 : PetscCall(VecSetType(ctx->v,bv->vtype));
452 268 : PetscCall(VecSetSizes(ctx->v,tlocal,PETSC_DECIDE));
453 268 : PetscCall(VecSetBlockSize(ctx->v,bs));
454 : }
455 348 : if (((PetscObject)bv)->name) {
456 122 : PetscCall(PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name));
457 122 : PetscCall(PetscObjectSetName((PetscObject)ctx->v,str));
458 : }
459 :
460 348 : if (PetscUnlikely(bv->Acreate)) {
461 2 : PetscCall(MatGetType(bv->Acreate,&mtype));
462 2 : PetscCall(PetscStrcmpAny(mtype,&isdense,MATSEQDENSE,MATMPIDENSE,""));
463 2 : PetscCheck(isdense,PetscObjectComm((PetscObject)bv->Acreate),PETSC_ERR_SUP,"BVSVEC requires a dense matrix in BVCreateFromMat()");
464 2 : PetscCall(MatDenseGetArrayRead(bv->Acreate,&aa));
465 2 : PetscCall(MatDenseGetLDA(bv->Acreate,&lda));
466 2 : PetscCall(VecGetArray(ctx->v,&vv));
467 18 : for (j=0;j<bv->m;j++) PetscCall(PetscArraycpy(vv+j*bv->ld,aa+j*lda,bv->n));
468 2 : PetscCall(VecRestoreArray(ctx->v,&vv));
469 2 : PetscCall(MatDenseRestoreArrayRead(bv->Acreate,&aa));
470 2 : PetscCall(MatDestroy(&bv->Acreate));
471 : }
472 :
473 348 : PetscCall(BVCreateVecEmpty(bv,&bv->cv[0]));
474 348 : PetscCall(BVCreateVecEmpty(bv,&bv->cv[1]));
475 :
476 348 : if (bv->cuda) {
477 : #if defined(PETSC_HAVE_CUDA)
478 : bv->ops->mult = BVMult_Svec_CUDA;
479 : bv->ops->multvec = BVMultVec_Svec_CUDA;
480 : bv->ops->multinplace = BVMultInPlace_Svec_CUDA;
481 : bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Svec_CUDA;
482 : bv->ops->dot = BVDot_Svec_CUDA;
483 : bv->ops->dotvec = BVDotVec_Svec_CUDA;
484 : bv->ops->dotvec_local = BVDotVec_Local_Svec_CUDA;
485 : bv->ops->scale = BVScale_Svec_CUDA;
486 : bv->ops->norm = BVNorm_Svec_CUDA;
487 : bv->ops->norm_local = BVNorm_Local_Svec_CUDA;
488 : bv->ops->normalize = BVNormalize_Svec_CUDA;
489 : bv->ops->matmult = BVMatMult_Svec_CUDA;
490 : bv->ops->copy = BVCopy_Svec_CUDA;
491 : bv->ops->copycolumn = BVCopyColumn_Svec_CUDA;
492 : bv->ops->resize = BVResize_Svec_CUDA;
493 : bv->ops->getcolumn = BVGetColumn_Svec_CUDA;
494 : bv->ops->restorecolumn = BVRestoreColumn_Svec_CUDA;
495 : bv->ops->restoresplit = BVRestoreSplit_Svec_CUDA;
496 : bv->ops->getmat = BVGetMat_Svec_CUDA;
497 : bv->ops->restoremat = BVRestoreMat_Svec_CUDA;
498 : #endif
499 : } else {
500 348 : bv->ops->mult = BVMult_Svec;
501 348 : bv->ops->multvec = BVMultVec_Svec;
502 348 : bv->ops->multinplace = BVMultInPlace_Svec;
503 348 : bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Svec;
504 348 : bv->ops->dot = BVDot_Svec;
505 348 : bv->ops->dotvec = BVDotVec_Svec;
506 348 : bv->ops->dotvec_local = BVDotVec_Local_Svec;
507 348 : bv->ops->scale = BVScale_Svec;
508 348 : bv->ops->norm = BVNorm_Svec;
509 348 : bv->ops->norm_local = BVNorm_Local_Svec;
510 348 : bv->ops->normalize = BVNormalize_Svec;
511 348 : bv->ops->matmult = BVMatMult_Svec;
512 348 : bv->ops->copy = BVCopy_Svec;
513 348 : bv->ops->copycolumn = BVCopyColumn_Svec;
514 348 : bv->ops->resize = BVResize_Svec;
515 348 : bv->ops->getcolumn = BVGetColumn_Svec;
516 348 : bv->ops->restorecolumn = BVRestoreColumn_Svec;
517 348 : bv->ops->getmat = BVGetMat_Default;
518 348 : bv->ops->restoremat = BVRestoreMat_Default;
519 : }
520 348 : bv->ops->getarray = BVGetArray_Svec;
521 348 : bv->ops->restorearray = BVRestoreArray_Svec;
522 348 : bv->ops->getarrayread = BVGetArrayRead_Svec;
523 348 : bv->ops->restorearrayread = BVRestoreArrayRead_Svec;
524 348 : bv->ops->destroy = BVDestroy_Svec;
525 348 : if (!ctx->mpi) bv->ops->view = BVView_Svec;
526 348 : PetscFunctionReturn(PETSC_SUCCESS);
527 : }
|