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 : SLEPc singular value solver: "cyclic"
12 :
13 : Method: Uses a Hermitian eigensolver for H(A) = [ 0 A ; A^T 0 ]
14 : */
15 :
16 : #include <slepc/private/svdimpl.h> /*I "slepcsvd.h" I*/
17 : #include <slepc/private/bvimpl.h>
18 : #include "cyclic.h"
19 :
20 2842 : static PetscErrorCode MatMult_Cyclic(Mat B,Vec x,Vec y)
21 : {
22 2842 : SVD_CYCLIC_SHELL *ctx;
23 2842 : const PetscScalar *px;
24 2842 : PetscScalar *py;
25 2842 : PetscInt m;
26 :
27 2842 : PetscFunctionBegin;
28 2842 : PetscCall(MatShellGetContext(B,&ctx));
29 2842 : PetscCall(MatGetLocalSize(ctx->A,&m,NULL));
30 2842 : PetscCall(VecGetArrayRead(x,&px));
31 2842 : PetscCall(VecGetArrayWrite(y,&py));
32 2842 : PetscCall(VecPlaceArray(ctx->x1,px));
33 2842 : PetscCall(VecPlaceArray(ctx->x2,px+m));
34 2842 : PetscCall(VecPlaceArray(ctx->y1,py));
35 2842 : PetscCall(VecPlaceArray(ctx->y2,py+m));
36 2842 : PetscCall(MatMult(ctx->A,ctx->x2,ctx->y1));
37 2842 : PetscCall(MatMult(ctx->AT,ctx->x1,ctx->y2));
38 2842 : PetscCall(VecResetArray(ctx->x1));
39 2842 : PetscCall(VecResetArray(ctx->x2));
40 2842 : PetscCall(VecResetArray(ctx->y1));
41 2842 : PetscCall(VecResetArray(ctx->y2));
42 2842 : PetscCall(VecRestoreArrayRead(x,&px));
43 2842 : PetscCall(VecRestoreArrayWrite(y,&py));
44 2842 : PetscFunctionReturn(PETSC_SUCCESS);
45 : }
46 :
47 0 : static PetscErrorCode MatGetDiagonal_Cyclic(Mat B,Vec diag)
48 : {
49 0 : PetscFunctionBegin;
50 0 : PetscCall(VecSet(diag,0.0));
51 0 : PetscFunctionReturn(PETSC_SUCCESS);
52 : }
53 :
54 27 : static PetscErrorCode MatDestroy_Cyclic(Mat B)
55 : {
56 27 : SVD_CYCLIC_SHELL *ctx;
57 :
58 27 : PetscFunctionBegin;
59 27 : PetscCall(MatShellGetContext(B,&ctx));
60 27 : PetscCall(VecDestroy(&ctx->x1));
61 27 : PetscCall(VecDestroy(&ctx->x2));
62 27 : PetscCall(VecDestroy(&ctx->y1));
63 27 : PetscCall(VecDestroy(&ctx->y2));
64 27 : if (ctx->misaligned) {
65 0 : PetscCall(VecDestroy(&ctx->wx2));
66 0 : PetscCall(VecDestroy(&ctx->wy2));
67 : }
68 27 : PetscCall(PetscFree(ctx));
69 27 : PetscFunctionReturn(PETSC_SUCCESS);
70 : }
71 :
72 : /*
73 : Builds cyclic matrix C = | 0 A |
74 : | AT 0 |
75 : */
76 42 : static PetscErrorCode SVDCyclicGetCyclicMat(SVD svd,Mat A,Mat AT,Mat *C)
77 : {
78 42 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
79 42 : SVD_CYCLIC_SHELL *ctx;
80 42 : PetscInt i,M,N,m,n,Istart,Iend;
81 42 : VecType vtype;
82 42 : Mat Zm,Zn;
83 : #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
84 : PetscBool gpu;
85 : const PetscInt *ranges;
86 : PetscMPIInt size;
87 : #endif
88 :
89 42 : PetscFunctionBegin;
90 42 : PetscCall(MatGetSize(A,&M,&N));
91 42 : PetscCall(MatGetLocalSize(A,&m,&n));
92 :
93 42 : if (cyclic->explicitmatrix) {
94 15 : PetscCheck(svd->expltrans,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Cannot use explicit cyclic matrix with implicit transpose");
95 15 : PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zm));
96 15 : PetscCall(MatSetSizes(Zm,m,m,M,M));
97 15 : PetscCall(MatSetFromOptions(Zm));
98 15 : PetscCall(MatGetOwnershipRange(Zm,&Istart,&Iend));
99 2957 : for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(Zm,i,i,0.0,INSERT_VALUES));
100 15 : PetscCall(MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY));
101 15 : PetscCall(MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY));
102 15 : PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zn));
103 15 : PetscCall(MatSetSizes(Zn,n,n,N,N));
104 15 : PetscCall(MatSetFromOptions(Zn));
105 15 : PetscCall(MatGetOwnershipRange(Zn,&Istart,&Iend));
106 1524 : for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(Zn,i,i,0.0,INSERT_VALUES));
107 15 : PetscCall(MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY));
108 15 : PetscCall(MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY));
109 15 : PetscCall(MatCreateTile(1.0,Zm,1.0,A,1.0,AT,1.0,Zn,C));
110 15 : PetscCall(MatDestroy(&Zm));
111 15 : PetscCall(MatDestroy(&Zn));
112 : } else {
113 27 : PetscCall(PetscNew(&ctx));
114 27 : ctx->A = A;
115 27 : ctx->AT = AT;
116 27 : ctx->swapped = svd->swapped;
117 27 : PetscCall(MatCreateVecsEmpty(A,&ctx->x2,&ctx->x1));
118 27 : PetscCall(MatCreateVecsEmpty(A,&ctx->y2,&ctx->y1));
119 27 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,ctx,C));
120 27 : PetscCall(MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cyclic));
121 27 : PetscCall(MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_Cyclic));
122 : #if defined(PETSC_HAVE_CUDA)
123 : PetscCall(PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&gpu,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,""));
124 : if (gpu) PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic_CUDA));
125 : else
126 : #elif defined(PETSC_HAVE_HIP)
127 : PetscCall(PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&gpu,MATSEQAIJHIPSPARSE,MATMPIAIJHIPSPARSE,""));
128 : if (gpu) PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic_HIP));
129 : else
130 : #endif
131 27 : PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic));
132 27 : PetscCall(MatGetVecType(A,&vtype));
133 27 : PetscCall(MatSetVecType(*C,vtype));
134 : #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
135 : if (gpu) {
136 : /* check alignment of bottom block */
137 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ctx->x1),&size));
138 : PetscCall(VecGetOwnershipRanges(ctx->x1,&ranges));
139 : for (i=0;i<size;i++) {
140 : ctx->misaligned = (((ranges[i+1]-ranges[i])*sizeof(PetscScalar))%16)? PETSC_TRUE: PETSC_FALSE;
141 : if (ctx->misaligned) break;
142 : }
143 : if (ctx->misaligned) { /* create work vectors for MatMult */
144 : PetscCall(VecDuplicate(ctx->x2,&ctx->wx2));
145 : PetscCall(VecDuplicate(ctx->y2,&ctx->wy2));
146 : }
147 : }
148 : #endif
149 : }
150 42 : PetscFunctionReturn(PETSC_SUCCESS);
151 : }
152 :
153 6635 : static PetscErrorCode MatMult_ECross(Mat B,Vec x,Vec y)
154 : {
155 6635 : SVD_CYCLIC_SHELL *ctx;
156 6635 : const PetscScalar *px;
157 6635 : PetscScalar *py;
158 6635 : PetscInt mn,m,n;
159 :
160 6635 : PetscFunctionBegin;
161 6635 : PetscCall(MatShellGetContext(B,&ctx));
162 6635 : PetscCall(MatGetLocalSize(ctx->A,NULL,&n));
163 6635 : PetscCall(VecGetLocalSize(y,&mn));
164 6635 : m = mn-n;
165 6635 : PetscCall(VecGetArrayRead(x,&px));
166 6635 : PetscCall(VecGetArrayWrite(y,&py));
167 6635 : PetscCall(VecPlaceArray(ctx->x1,px));
168 6635 : PetscCall(VecPlaceArray(ctx->x2,px+m));
169 6635 : PetscCall(VecPlaceArray(ctx->y1,py));
170 6635 : PetscCall(VecPlaceArray(ctx->y2,py+m));
171 6635 : PetscCall(VecCopy(ctx->x1,ctx->y1));
172 6635 : PetscCall(MatMult(ctx->A,ctx->x2,ctx->w));
173 6635 : PetscCall(MatMult(ctx->AT,ctx->w,ctx->y2));
174 6635 : PetscCall(VecResetArray(ctx->x1));
175 6635 : PetscCall(VecResetArray(ctx->x2));
176 6635 : PetscCall(VecResetArray(ctx->y1));
177 6635 : PetscCall(VecResetArray(ctx->y2));
178 6635 : PetscCall(VecRestoreArrayRead(x,&px));
179 6635 : PetscCall(VecRestoreArrayWrite(y,&py));
180 6635 : PetscFunctionReturn(PETSC_SUCCESS);
181 : }
182 :
183 5 : static PetscErrorCode MatGetDiagonal_ECross(Mat B,Vec d)
184 : {
185 5 : SVD_CYCLIC_SHELL *ctx;
186 5 : PetscScalar *pd;
187 5 : PetscMPIInt len;
188 5 : PetscInt mn,m,n,N,i,j,start,end,ncols;
189 5 : PetscScalar *work1,*work2,*diag;
190 5 : const PetscInt *cols;
191 5 : const PetscScalar *vals;
192 :
193 5 : PetscFunctionBegin;
194 5 : PetscCall(MatShellGetContext(B,&ctx));
195 5 : PetscCall(MatGetLocalSize(ctx->A,NULL,&n));
196 5 : PetscCall(VecGetLocalSize(d,&mn));
197 5 : m = mn-n;
198 5 : PetscCall(VecGetArrayWrite(d,&pd));
199 5 : PetscCall(VecPlaceArray(ctx->y1,pd));
200 5 : PetscCall(VecSet(ctx->y1,1.0));
201 5 : PetscCall(VecResetArray(ctx->y1));
202 5 : PetscCall(VecPlaceArray(ctx->y2,pd+m));
203 5 : if (!ctx->diag) {
204 : /* compute diagonal from rows and store in ctx->diag */
205 5 : PetscCall(VecDuplicate(ctx->y2,&ctx->diag));
206 5 : PetscCall(MatGetSize(ctx->A,NULL,&N));
207 5 : PetscCall(PetscCalloc2(N,&work1,N,&work2));
208 5 : if (ctx->swapped) {
209 0 : PetscCall(MatGetOwnershipRange(ctx->AT,&start,&end));
210 0 : for (i=start;i<end;i++) {
211 0 : PetscCall(MatGetRow(ctx->AT,i,&ncols,NULL,&vals));
212 0 : for (j=0;j<ncols;j++) work1[i] += vals[j]*vals[j];
213 0 : PetscCall(MatRestoreRow(ctx->AT,i,&ncols,NULL,&vals));
214 : }
215 : } else {
216 5 : PetscCall(MatGetOwnershipRange(ctx->A,&start,&end));
217 111 : for (i=start;i<end;i++) {
218 106 : PetscCall(MatGetRow(ctx->A,i,&ncols,&cols,&vals));
219 806 : for (j=0;j<ncols;j++) work1[cols[j]] += vals[j]*vals[j];
220 106 : PetscCall(MatRestoreRow(ctx->A,i,&ncols,&cols,&vals));
221 : }
222 : }
223 5 : PetscCall(PetscMPIIntCast(N,&len));
224 15 : PetscCallMPI(MPIU_Allreduce(work1,work2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)B)));
225 5 : PetscCall(VecGetOwnershipRange(ctx->diag,&start,&end));
226 5 : PetscCall(VecGetArrayWrite(ctx->diag,&diag));
227 100 : for (i=start;i<end;i++) diag[i-start] = work2[i];
228 5 : PetscCall(VecRestoreArrayWrite(ctx->diag,&diag));
229 5 : PetscCall(PetscFree2(work1,work2));
230 : }
231 5 : PetscCall(VecCopy(ctx->diag,ctx->y2));
232 5 : PetscCall(VecResetArray(ctx->y2));
233 5 : PetscCall(VecRestoreArrayWrite(d,&pd));
234 5 : PetscFunctionReturn(PETSC_SUCCESS);
235 : }
236 :
237 5 : static PetscErrorCode MatDestroy_ECross(Mat B)
238 : {
239 5 : SVD_CYCLIC_SHELL *ctx;
240 :
241 5 : PetscFunctionBegin;
242 5 : PetscCall(MatShellGetContext(B,&ctx));
243 5 : PetscCall(VecDestroy(&ctx->x1));
244 5 : PetscCall(VecDestroy(&ctx->x2));
245 5 : PetscCall(VecDestroy(&ctx->y1));
246 5 : PetscCall(VecDestroy(&ctx->y2));
247 5 : PetscCall(VecDestroy(&ctx->diag));
248 5 : PetscCall(VecDestroy(&ctx->w));
249 5 : if (ctx->misaligned) {
250 0 : PetscCall(VecDestroy(&ctx->wx2));
251 0 : PetscCall(VecDestroy(&ctx->wy2));
252 : }
253 5 : PetscCall(PetscFree(ctx));
254 5 : PetscFunctionReturn(PETSC_SUCCESS);
255 : }
256 :
257 : /*
258 : Builds extended cross product matrix C = | I_m 0 |
259 : | 0 AT*A |
260 : t is an auxiliary Vec used to take the dimensions of the upper block
261 : */
262 10 : static PetscErrorCode SVDCyclicGetECrossMat(SVD svd,Mat A,Mat AT,Mat *C,Vec t)
263 : {
264 10 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
265 10 : SVD_CYCLIC_SHELL *ctx;
266 10 : PetscInt i,M,N,m,n,Istart,Iend;
267 10 : VecType vtype;
268 10 : Mat Id,Zm,Zn,ATA;
269 : #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
270 : PetscBool gpu;
271 : const PetscInt *ranges;
272 : PetscMPIInt size;
273 : #endif
274 :
275 10 : PetscFunctionBegin;
276 10 : PetscCall(MatGetSize(A,NULL,&N));
277 10 : PetscCall(MatGetLocalSize(A,NULL,&n));
278 10 : PetscCall(VecGetSize(t,&M));
279 10 : PetscCall(VecGetLocalSize(t,&m));
280 :
281 10 : if (cyclic->explicitmatrix) {
282 5 : PetscCheck(svd->expltrans,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Cannot use explicit cyclic matrix with implicit transpose");
283 5 : PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)svd),m,m,M,M,1.0,&Id));
284 5 : PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zm));
285 5 : PetscCall(MatSetSizes(Zm,m,n,M,N));
286 5 : PetscCall(MatSetFromOptions(Zm));
287 5 : PetscCall(MatGetOwnershipRange(Zm,&Istart,&Iend));
288 394 : for (i=Istart;i<Iend;i++) {
289 389 : if (i<N) PetscCall(MatSetValue(Zm,i,i,0.0,INSERT_VALUES));
290 : }
291 5 : PetscCall(MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY));
292 5 : PetscCall(MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY));
293 5 : PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zn));
294 5 : PetscCall(MatSetSizes(Zn,n,m,N,M));
295 5 : PetscCall(MatSetFromOptions(Zn));
296 5 : PetscCall(MatGetOwnershipRange(Zn,&Istart,&Iend));
297 404 : for (i=Istart;i<Iend;i++) {
298 399 : if (i<m) PetscCall(MatSetValue(Zn,i,i,0.0,INSERT_VALUES));
299 : }
300 5 : PetscCall(MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY));
301 5 : PetscCall(MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY));
302 5 : PetscCall(MatProductCreate(AT,A,NULL,&ATA));
303 5 : PetscCall(MatProductSetType(ATA,MATPRODUCT_AB));
304 5 : PetscCall(MatProductSetFromOptions(ATA));
305 5 : PetscCall(MatProductSymbolic(ATA));
306 5 : PetscCall(MatProductNumeric(ATA));
307 5 : PetscCall(MatCreateTile(1.0,Id,1.0,Zm,1.0,Zn,1.0,ATA,C));
308 5 : PetscCall(MatDestroy(&Id));
309 5 : PetscCall(MatDestroy(&Zm));
310 5 : PetscCall(MatDestroy(&Zn));
311 5 : PetscCall(MatDestroy(&ATA));
312 : } else {
313 5 : PetscCall(PetscNew(&ctx));
314 5 : ctx->A = A;
315 5 : ctx->AT = AT;
316 5 : ctx->swapped = svd->swapped;
317 5 : PetscCall(VecDuplicateEmpty(t,&ctx->x1));
318 5 : PetscCall(VecDuplicateEmpty(t,&ctx->y1));
319 5 : PetscCall(MatCreateVecsEmpty(A,&ctx->x2,NULL));
320 5 : PetscCall(MatCreateVecsEmpty(A,&ctx->y2,NULL));
321 5 : PetscCall(MatCreateVecs(A,NULL,&ctx->w));
322 5 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,ctx,C));
323 5 : PetscCall(MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ECross));
324 5 : PetscCall(MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_ECross));
325 : #if defined(PETSC_HAVE_CUDA)
326 : PetscCall(PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&gpu,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,""));
327 : if (gpu) PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross_CUDA));
328 : else
329 : #elif defined(PETSC_HAVE_HIP)
330 : PetscCall(PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&gpu,MATSEQAIJHIPSPARSE,MATMPIAIJHIPSPARSE,""));
331 : if (gpu) PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross_HIP));
332 : else
333 : #endif
334 5 : PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross));
335 5 : PetscCall(MatGetVecType(A,&vtype));
336 5 : PetscCall(MatSetVecType(*C,vtype));
337 : #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
338 : if (gpu) {
339 : /* check alignment of bottom block */
340 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ctx->x1),&size));
341 : PetscCall(VecGetOwnershipRanges(ctx->x1,&ranges));
342 : for (i=0;i<size;i++) {
343 : ctx->misaligned = (((ranges[i+1]-ranges[i])*sizeof(PetscScalar))%16)? PETSC_TRUE: PETSC_FALSE;
344 : if (ctx->misaligned) break;
345 : }
346 : if (ctx->misaligned) { /* create work vectors for MatMult */
347 : PetscCall(VecDuplicate(ctx->x2,&ctx->wx2));
348 : PetscCall(VecDuplicate(ctx->y2,&ctx->wy2));
349 : }
350 : }
351 : #endif
352 : }
353 10 : PetscFunctionReturn(PETSC_SUCCESS);
354 : }
355 :
356 : /* Convergence test relative to the norm of R (used in GSVD only) */
357 136 : static PetscErrorCode EPSConv_Cyclic(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
358 : {
359 136 : SVD svd = (SVD)ctx;
360 :
361 136 : PetscFunctionBegin;
362 136 : *errest = res/PetscMax(svd->nrma,svd->nrmb);
363 136 : PetscFunctionReturn(PETSC_SUCCESS);
364 : }
365 :
366 42 : static PetscErrorCode SVDSetUp_Cyclic(SVD svd)
367 : {
368 42 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
369 42 : PetscInt M,N,m,n,p,k,i,isl,offset,nev,ncv,mpd,maxit;
370 42 : PetscReal tol;
371 42 : const PetscScalar *isa,*oa;
372 42 : PetscScalar *va;
373 42 : EPSProblemType ptype;
374 42 : PetscBool trackall,issinv;
375 42 : Vec v,t;
376 42 : ST st;
377 42 : Mat Omega;
378 42 : MatType Atype;
379 :
380 42 : PetscFunctionBegin;
381 42 : if (svd->nsv==0) svd->nsv = 1;
382 42 : PetscCall(MatGetSize(svd->A,&M,&N));
383 42 : PetscCall(MatGetLocalSize(svd->A,&m,&n));
384 42 : if (!cyclic->eps) PetscCall(SVDCyclicGetEPS(svd,&cyclic->eps));
385 42 : PetscCall(MatDestroy(&cyclic->C));
386 42 : PetscCall(MatDestroy(&cyclic->D));
387 42 : if (svd->isgeneralized) {
388 10 : if (svd->which==SVD_SMALLEST) { /* alternative pencil */
389 0 : PetscCall(MatCreateVecs(svd->B,NULL,&t));
390 0 : PetscCall(SVDCyclicGetCyclicMat(svd,svd->B,svd->BT,&cyclic->C));
391 0 : PetscCall(SVDCyclicGetECrossMat(svd,svd->A,svd->AT,&cyclic->D,t));
392 : } else {
393 10 : PetscCall(MatCreateVecs(svd->A,NULL,&t));
394 10 : PetscCall(SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C));
395 10 : PetscCall(SVDCyclicGetECrossMat(svd,svd->B,svd->BT,&cyclic->D,t));
396 : }
397 10 : PetscCall(VecDestroy(&t));
398 10 : PetscCall(EPSSetOperators(cyclic->eps,cyclic->C,cyclic->D));
399 10 : PetscCall(EPSGetProblemType(cyclic->eps,&ptype));
400 10 : if (!ptype) PetscCall(EPSSetProblemType(cyclic->eps,EPS_GHEP));
401 32 : } else if (svd->ishyperbolic) {
402 7 : PetscCall(SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C));
403 7 : PetscCall(MatCreateVecs(cyclic->C,&v,NULL));
404 7 : PetscCall(VecSet(v,1.0));
405 7 : PetscCall(VecGetArrayRead(svd->omega,&oa));
406 7 : PetscCall(VecGetArray(v,&va));
407 7 : if (svd->swapped) PetscCall(PetscArraycpy(va+m,oa,n));
408 5 : else PetscCall(PetscArraycpy(va,oa,m));
409 7 : PetscCall(VecRestoreArrayRead(svd->omega,&oa));
410 7 : PetscCall(VecRestoreArray(v,&va));
411 7 : PetscCall(MatGetType(svd->OP,&Atype));
412 7 : PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Omega));
413 7 : PetscCall(MatSetSizes(Omega,m+n,m+n,M+N,M+N));
414 7 : PetscCall(MatSetType(Omega,Atype));
415 7 : PetscCall(MatDiagonalSet(Omega,v,INSERT_VALUES));
416 7 : PetscCall(EPSSetOperators(cyclic->eps,cyclic->C,Omega));
417 7 : PetscCall(EPSSetProblemType(cyclic->eps,EPS_GHIEP));
418 7 : PetscCall(MatDestroy(&Omega));
419 7 : PetscCall(VecDestroy(&v));
420 : } else {
421 25 : PetscCall(SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C));
422 25 : PetscCall(EPSSetOperators(cyclic->eps,cyclic->C,NULL));
423 25 : PetscCall(EPSSetProblemType(cyclic->eps,EPS_HEP));
424 : }
425 42 : if (!cyclic->usereps) {
426 41 : if (svd->which == SVD_LARGEST) {
427 39 : PetscCall(EPSGetST(cyclic->eps,&st));
428 39 : PetscCall(PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv));
429 39 : if (issinv) PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_TARGET_MAGNITUDE));
430 38 : else if (svd->ishyperbolic) PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_MAGNITUDE));
431 31 : else PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL));
432 : } else {
433 2 : if (svd->isgeneralized) { /* computes sigma^{-1} via alternative pencil */
434 0 : PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL));
435 : } else {
436 2 : if (svd->ishyperbolic) PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_TARGET_MAGNITUDE));
437 2 : else PetscCall(EPSSetEigenvalueComparison(cyclic->eps,SlepcCompareSmallestPosReal,NULL));
438 2 : PetscCall(EPSSetTarget(cyclic->eps,0.0));
439 : }
440 : }
441 41 : PetscCall(EPSGetDimensions(cyclic->eps,&nev,&ncv,&mpd));
442 41 : PetscCheck(nev==1 || nev>=2*svd->nsv,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"The number of requested eigenvalues %" PetscInt_FMT " must be at least 2*%" PetscInt_FMT,nev,svd->nsv);
443 41 : nev = PetscMax(nev,2*svd->nsv);
444 41 : if (ncv==PETSC_DETERMINE && svd->ncv!=PETSC_DETERMINE) ncv = PetscMax(3*svd->nsv,svd->ncv);
445 41 : if (mpd==PETSC_DETERMINE && svd->mpd!=PETSC_DETERMINE) mpd = svd->mpd;
446 41 : PetscCall(EPSSetDimensions(cyclic->eps,nev,ncv,mpd));
447 41 : PetscCall(EPSGetTolerances(cyclic->eps,&tol,&maxit));
448 47 : if (tol==(PetscReal)PETSC_DETERMINE) tol = svd->tol==(PetscReal)PETSC_DETERMINE? SLEPC_DEFAULT_TOL/10.0: svd->tol;
449 41 : if (maxit==PETSC_DETERMINE) maxit = svd->max_it;
450 41 : PetscCall(EPSSetTolerances(cyclic->eps,tol,maxit));
451 41 : switch (svd->conv) {
452 3 : case SVD_CONV_ABS:
453 3 : PetscCall(EPSSetConvergenceTest(cyclic->eps,EPS_CONV_ABS));break;
454 28 : case SVD_CONV_REL:
455 28 : PetscCall(EPSSetConvergenceTest(cyclic->eps,EPS_CONV_REL));break;
456 10 : case SVD_CONV_NORM:
457 10 : if (svd->isgeneralized) {
458 10 : if (!svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
459 10 : if (!svd->nrmb) PetscCall(MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb));
460 10 : PetscCall(EPSSetConvergenceTestFunction(cyclic->eps,EPSConv_Cyclic,svd,NULL));
461 : } else {
462 0 : PetscCall(EPSSetConvergenceTest(cyclic->eps,EPS_CONV_NORM));break;
463 : }
464 : break;
465 0 : case SVD_CONV_MAXIT:
466 0 : SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Maxit convergence test not supported in this solver");
467 0 : case SVD_CONV_USER:
468 0 : SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
469 : }
470 1 : }
471 42 : SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
472 : /* Transfer the trackall option from svd to eps */
473 42 : PetscCall(SVDGetTrackAll(svd,&trackall));
474 42 : PetscCall(EPSSetTrackAll(cyclic->eps,trackall));
475 : /* Transfer the initial subspace from svd to eps */
476 42 : if (svd->nini<0 || svd->ninil<0) {
477 8 : for (i=0;i<-PetscMin(svd->nini,svd->ninil);i++) {
478 4 : PetscCall(MatCreateVecs(cyclic->C,&v,NULL));
479 4 : PetscCall(VecGetArrayWrite(v,&va));
480 4 : if (svd->isgeneralized) PetscCall(MatGetLocalSize(svd->B,&p,NULL));
481 4 : k = (svd->isgeneralized && svd->which==SVD_SMALLEST)? p: m; /* size of upper block row */
482 4 : if (i<-svd->ninil) {
483 4 : PetscCall(VecGetArrayRead(svd->ISL[i],&isa));
484 4 : if (svd->isgeneralized) {
485 1 : PetscCall(VecGetLocalSize(svd->ISL[i],&isl));
486 1 : PetscCheck(isl==m+p,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Size mismatch for left initial vector");
487 1 : offset = (svd->which==SVD_SMALLEST)? m: 0;
488 1 : PetscCall(PetscArraycpy(va,isa+offset,k));
489 : } else {
490 3 : PetscCall(VecGetLocalSize(svd->ISL[i],&isl));
491 3 : PetscCheck(isl==k,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Size mismatch for left initial vector");
492 3 : PetscCall(PetscArraycpy(va,isa,k));
493 : }
494 4 : PetscCall(VecRestoreArrayRead(svd->IS[i],&isa));
495 0 : } else PetscCall(PetscArrayzero(&va,k));
496 4 : if (i<-svd->nini) {
497 4 : PetscCall(VecGetLocalSize(svd->IS[i],&isl));
498 4 : PetscCheck(isl==n,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Size mismatch for right initial vector");
499 4 : PetscCall(VecGetArrayRead(svd->IS[i],&isa));
500 4 : PetscCall(PetscArraycpy(va+k,isa,n));
501 4 : PetscCall(VecRestoreArrayRead(svd->IS[i],&isa));
502 0 : } else PetscCall(PetscArrayzero(va+k,n));
503 4 : PetscCall(VecRestoreArrayWrite(v,&va));
504 4 : PetscCall(VecDestroy(&svd->IS[i]));
505 4 : svd->IS[i] = v;
506 : }
507 4 : svd->nini = PetscMin(svd->nini,svd->ninil);
508 4 : PetscCall(EPSSetInitialSpace(cyclic->eps,-svd->nini,svd->IS));
509 4 : PetscCall(SlepcBasisDestroy_Private(&svd->nini,&svd->IS));
510 4 : PetscCall(SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL));
511 : }
512 42 : PetscCall(EPSSetUp(cyclic->eps));
513 42 : PetscCall(EPSGetDimensions(cyclic->eps,NULL,&svd->ncv,&svd->mpd));
514 42 : svd->ncv = PetscMin(svd->ncv,PetscMin(M,N));
515 42 : PetscCall(EPSGetTolerances(cyclic->eps,NULL,&svd->max_it));
516 42 : if (svd->tol==(PetscReal)PETSC_DETERMINE) svd->tol = SLEPC_DEFAULT_TOL;
517 :
518 42 : svd->leftbasis = PETSC_TRUE;
519 42 : PetscCall(SVDAllocateSolution(svd,0));
520 42 : PetscFunctionReturn(PETSC_SUCCESS);
521 : }
522 :
523 8497 : static PetscErrorCode SVDCyclicCheckEigenvalue(SVD svd,PetscScalar er,PetscScalar ei,PetscReal *sigma,PetscBool *isreal)
524 : {
525 8497 : PetscFunctionBegin;
526 8497 : if (svd->ishyperbolic && PetscDefined(USE_COMPLEX) && PetscAbsReal(PetscImaginaryPart(er))>10*PetscAbsReal(PetscRealPart(er))) {
527 228 : *sigma = PetscImaginaryPart(er);
528 228 : if (isreal) *isreal = PETSC_FALSE;
529 8269 : } else if (svd->ishyperbolic && !PetscDefined(USE_COMPLEX) && PetscAbsScalar(ei)>10*PetscAbsScalar(er)) {
530 : *sigma = PetscRealPart(ei);
531 : if (isreal) *isreal = PETSC_FALSE;
532 : } else {
533 8269 : *sigma = PetscRealPart(er);
534 8269 : if (isreal) *isreal = PETSC_TRUE;
535 : }
536 8497 : PetscFunctionReturn(PETSC_SUCCESS);
537 : }
538 :
539 42 : static PetscErrorCode SVDSolve_Cyclic(SVD svd)
540 : {
541 42 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
542 42 : PetscInt i,j,nconv;
543 42 : PetscScalar er,ei;
544 42 : PetscReal sigma;
545 :
546 42 : PetscFunctionBegin;
547 42 : PetscCall(EPSSolve(cyclic->eps));
548 42 : PetscCall(EPSGetConverged(cyclic->eps,&nconv));
549 42 : PetscCall(EPSGetIterationNumber(cyclic->eps,&svd->its));
550 42 : PetscCall(EPSGetConvergedReason(cyclic->eps,(EPSConvergedReason*)&svd->reason));
551 671 : for (i=0,j=0;i<nconv;i++) {
552 629 : PetscCall(EPSGetEigenvalue(cyclic->eps,i,&er,&ei));
553 629 : PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL));
554 629 : if (sigma>0.0) {
555 428 : if (svd->isgeneralized && svd->which==SVD_SMALLEST) svd->sigma[j] = 1.0/sigma;
556 428 : else svd->sigma[j] = sigma;
557 428 : j++;
558 : }
559 : }
560 42 : svd->nconv = j;
561 42 : PetscFunctionReturn(PETSC_SUCCESS);
562 : }
563 :
564 21 : static PetscErrorCode SVDComputeVectors_Cyclic_Standard(SVD svd)
565 : {
566 21 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
567 21 : PetscInt i,j,m,nconv;
568 21 : PetscScalar er,ei;
569 21 : PetscReal sigma;
570 21 : const PetscScalar *px;
571 21 : Vec x,x1,x2;
572 :
573 21 : PetscFunctionBegin;
574 21 : PetscCall(MatCreateVecs(cyclic->C,&x,NULL));
575 21 : PetscCall(MatGetLocalSize(svd->A,&m,NULL));
576 21 : PetscCall(MatCreateVecsEmpty(svd->A,&x2,&x1));
577 21 : PetscCall(EPSGetConverged(cyclic->eps,&nconv));
578 188 : for (i=0,j=0;i<nconv;i++) {
579 167 : PetscCall(EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,NULL));
580 167 : PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL));
581 167 : if (sigma<0.0) continue;
582 155 : PetscCall(VecGetArrayRead(x,&px));
583 155 : PetscCall(VecPlaceArray(x1,px));
584 155 : PetscCall(VecPlaceArray(x2,px+m));
585 155 : PetscCall(BVInsertVec(svd->U,j,x1));
586 155 : PetscCall(BVScaleColumn(svd->U,j,PETSC_SQRT2));
587 155 : PetscCall(BVInsertVec(svd->V,j,x2));
588 155 : PetscCall(BVScaleColumn(svd->V,j,PETSC_SQRT2));
589 155 : PetscCall(VecResetArray(x1));
590 155 : PetscCall(VecResetArray(x2));
591 155 : PetscCall(VecRestoreArrayRead(x,&px));
592 155 : j++;
593 : }
594 21 : PetscCall(VecDestroy(&x));
595 21 : PetscCall(VecDestroy(&x1));
596 21 : PetscCall(VecDestroy(&x2));
597 21 : PetscFunctionReturn(PETSC_SUCCESS);
598 : }
599 :
600 10 : static PetscErrorCode SVDComputeVectors_Cyclic_Generalized(SVD svd)
601 : {
602 10 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
603 10 : PetscInt i,j,m,p,nconv;
604 10 : PetscScalar *dst,er,ei;
605 10 : PetscReal sigma;
606 10 : const PetscScalar *src,*px;
607 10 : Vec u,v,x,x1,x2,uv;
608 :
609 10 : PetscFunctionBegin;
610 10 : PetscCall(MatGetLocalSize(svd->A,&m,NULL));
611 10 : PetscCall(MatGetLocalSize(svd->B,&p,NULL));
612 10 : PetscCall(MatCreateVecs(cyclic->C,&x,NULL));
613 10 : if (svd->which==SVD_SMALLEST) PetscCall(MatCreateVecsEmpty(svd->B,&x1,&x2));
614 10 : else PetscCall(MatCreateVecsEmpty(svd->A,&x2,&x1));
615 10 : PetscCall(MatCreateVecs(svd->A,NULL,&u));
616 10 : PetscCall(MatCreateVecs(svd->B,NULL,&v));
617 10 : PetscCall(EPSGetConverged(cyclic->eps,&nconv));
618 80 : for (i=0,j=0;i<nconv;i++) {
619 70 : PetscCall(EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,NULL));
620 70 : PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL));
621 70 : if (sigma<0.0) continue;
622 70 : if (svd->which==SVD_SMALLEST) {
623 : /* evec_i = 1/sqrt(2)*[ v_i; w_i ], w_i = x_i/c_i */
624 0 : PetscCall(VecGetArrayRead(x,&px));
625 0 : PetscCall(VecPlaceArray(x2,px));
626 0 : PetscCall(VecPlaceArray(x1,px+p));
627 0 : PetscCall(VecCopy(x2,v));
628 0 : PetscCall(VecScale(v,PETSC_SQRT2)); /* v_i = sqrt(2)*evec_i_1 */
629 0 : PetscCall(VecScale(x1,PETSC_SQRT2)); /* w_i = sqrt(2)*evec_i_2 */
630 0 : PetscCall(MatMult(svd->A,x1,u)); /* A*w_i = u_i */
631 0 : PetscCall(VecScale(x1,1.0/PetscSqrtScalar(1.0+sigma*sigma))); /* x_i = w_i*c_i */
632 0 : PetscCall(BVInsertVec(svd->V,j,x1));
633 0 : PetscCall(VecResetArray(x2));
634 0 : PetscCall(VecResetArray(x1));
635 0 : PetscCall(VecRestoreArrayRead(x,&px));
636 : } else {
637 : /* evec_i = 1/sqrt(2)*[ u_i; w_i ], w_i = x_i/s_i */
638 70 : PetscCall(VecGetArrayRead(x,&px));
639 70 : PetscCall(VecPlaceArray(x1,px));
640 70 : PetscCall(VecPlaceArray(x2,px+m));
641 70 : PetscCall(VecCopy(x1,u));
642 70 : PetscCall(VecScale(u,PETSC_SQRT2)); /* u_i = sqrt(2)*evec_i_1 */
643 70 : PetscCall(VecScale(x2,PETSC_SQRT2)); /* w_i = sqrt(2)*evec_i_2 */
644 70 : PetscCall(MatMult(svd->B,x2,v)); /* B*w_i = v_i */
645 70 : PetscCall(VecScale(x2,1.0/PetscSqrtScalar(1.0+sigma*sigma))); /* x_i = w_i*s_i */
646 70 : PetscCall(BVInsertVec(svd->V,j,x2));
647 70 : PetscCall(VecResetArray(x1));
648 70 : PetscCall(VecResetArray(x2));
649 70 : PetscCall(VecRestoreArrayRead(x,&px));
650 : }
651 : /* copy [u;v] to U[j] */
652 70 : PetscCall(BVGetColumn(svd->U,j,&uv));
653 70 : PetscCall(VecGetArrayWrite(uv,&dst));
654 70 : PetscCall(VecGetArrayRead(u,&src));
655 70 : PetscCall(PetscArraycpy(dst,src,m));
656 70 : PetscCall(VecRestoreArrayRead(u,&src));
657 70 : PetscCall(VecGetArrayRead(v,&src));
658 70 : PetscCall(PetscArraycpy(dst+m,src,p));
659 70 : PetscCall(VecRestoreArrayRead(v,&src));
660 70 : PetscCall(VecRestoreArrayWrite(uv,&dst));
661 70 : PetscCall(BVRestoreColumn(svd->U,j,&uv));
662 70 : j++;
663 : }
664 10 : PetscCall(VecDestroy(&x));
665 10 : PetscCall(VecDestroy(&x1));
666 10 : PetscCall(VecDestroy(&x2));
667 10 : PetscCall(VecDestroy(&u));
668 10 : PetscCall(VecDestroy(&v));
669 10 : PetscFunctionReturn(PETSC_SUCCESS);
670 : }
671 :
672 : #if defined(PETSC_USE_COMPLEX)
673 : /* VecMaxAbs: returns the entry of x that has max(abs(x(i))), using w as a workspace vector */
674 4 : static PetscErrorCode VecMaxAbs(Vec x,Vec w,PetscScalar *v)
675 : {
676 4 : PetscMPIInt size,rank,root;
677 4 : const PetscScalar *xx;
678 4 : const PetscInt *ranges;
679 4 : PetscReal val;
680 4 : PetscInt p;
681 :
682 4 : PetscFunctionBegin;
683 4 : PetscCall(VecCopy(x,w));
684 4 : PetscCall(VecAbs(w));
685 4 : PetscCall(VecMax(w,&p,&val));
686 4 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)x),&size));
687 4 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)x),&rank));
688 4 : PetscCall(VecGetOwnershipRanges(x,&ranges));
689 4 : for (root=0;root<size;root++) if (p>=ranges[root] && p<ranges[root+1]) break;
690 4 : if (rank==root) {
691 4 : PetscCall(VecGetArrayRead(x,&xx));
692 4 : *v = xx[p-ranges[root]];
693 4 : PetscCall(VecRestoreArrayRead(x,&xx));
694 : }
695 8 : PetscCallMPI(MPI_Bcast(v,1,MPIU_SCALAR,root,PetscObjectComm((PetscObject)x)));
696 4 : PetscFunctionReturn(PETSC_SUCCESS);
697 : }
698 : #endif
699 :
700 7 : static PetscErrorCode SVDComputeVectors_Cyclic_Hyperbolic(SVD svd)
701 : {
702 7 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
703 7 : PetscInt i,j,m,n,nconv;
704 7 : PetscScalar er,ei;
705 7 : PetscReal sigma,nrm;
706 7 : PetscBool isreal;
707 7 : const PetscScalar *px;
708 7 : Vec u,x,xi=NULL,x1,x2,x1i=NULL,x2i;
709 7 : BV U=NULL,V=NULL;
710 : #if !defined(PETSC_USE_COMPLEX)
711 : const PetscScalar *pxi;
712 : PetscReal nrmr,nrmi;
713 : #else
714 7 : PetscScalar alpha;
715 : #endif
716 :
717 7 : PetscFunctionBegin;
718 7 : PetscCall(MatCreateVecs(cyclic->C,&x,svd->ishyperbolic?&xi:NULL));
719 7 : PetscCall(MatGetLocalSize(svd->A,&m,NULL));
720 7 : PetscCall(MatCreateVecsEmpty(svd->OP,&x2,&x1));
721 : #if defined(PETSC_USE_COMPLEX)
722 7 : PetscCall(MatCreateVecs(svd->OP,&x2i,&x1i));
723 : #else
724 : PetscCall(MatCreateVecsEmpty(svd->OP,&x2i,&x1i));
725 : #endif
726 : /* set-up Omega-normalization of U */
727 7 : U = svd->swapped? svd->V: svd->U;
728 7 : V = svd->swapped? svd->U: svd->V;
729 7 : PetscCall(BVGetSizes(U,&n,NULL,NULL));
730 7 : PetscCall(BV_SetMatrixDiagonal(U,svd->omega,svd->A));
731 7 : PetscCall(EPSGetConverged(cyclic->eps,&nconv));
732 388 : for (i=0,j=0;i<nconv;i++) {
733 381 : PetscCall(EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,xi));
734 381 : PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,&isreal));
735 381 : if (sigma<0.0) continue;
736 192 : PetscCall(VecGetArrayRead(x,&px));
737 192 : if (svd->swapped) {
738 14 : PetscCall(VecPlaceArray(x2,px));
739 14 : PetscCall(VecPlaceArray(x1,px+m));
740 : } else {
741 178 : PetscCall(VecPlaceArray(x1,px));
742 178 : PetscCall(VecPlaceArray(x2,px+n));
743 : }
744 : #if defined(PETSC_USE_COMPLEX)
745 192 : PetscCall(BVInsertVec(U,j,x1));
746 192 : PetscCall(BVInsertVec(V,j,x2));
747 192 : if (!isreal) {
748 4 : PetscCall(VecMaxAbs(x1,x1i,&alpha));
749 4 : PetscCall(BVScaleColumn(U,j,PetscAbsScalar(alpha)/alpha));
750 4 : PetscCall(BVScaleColumn(V,j,PetscAbsScalar(alpha)/(alpha*PETSC_i)));
751 : }
752 : #else
753 : PetscCall(VecGetArrayRead(xi,&pxi));
754 : if (svd->swapped) {
755 : PetscCall(VecPlaceArray(x2i,pxi));
756 : PetscCall(VecPlaceArray(x1i,pxi+m));
757 : } else {
758 : PetscCall(VecPlaceArray(x1i,pxi));
759 : PetscCall(VecPlaceArray(x2i,pxi+n));
760 : }
761 : PetscCall(VecNorm(x2,NORM_2,&nrmr));
762 : PetscCall(VecNorm(x2i,NORM_2,&nrmi));
763 : if (nrmi>nrmr) {
764 : if (isreal) {
765 : PetscCall(BVInsertVec(U,j,x1i));
766 : PetscCall(BVInsertVec(V,j,x2i));
767 : } else {
768 : PetscCall(BVInsertVec(U,j,x1));
769 : PetscCall(BVInsertVec(V,j,x2i));
770 : }
771 : } else {
772 : if (isreal) {
773 : PetscCall(BVInsertVec(U,j,x1));
774 : PetscCall(BVInsertVec(V,j,x2));
775 : } else {
776 : PetscCall(BVInsertVec(U,j,x1i));
777 : PetscCall(BVScaleColumn(U,j,-1.0));
778 : PetscCall(BVInsertVec(V,j,x2));
779 : }
780 : }
781 : PetscCall(VecResetArray(x1i));
782 : PetscCall(VecResetArray(x2i));
783 : PetscCall(VecRestoreArrayRead(xi,&pxi));
784 : #endif
785 192 : PetscCall(VecResetArray(x1));
786 192 : PetscCall(VecResetArray(x2));
787 192 : PetscCall(VecRestoreArrayRead(x,&px));
788 192 : PetscCall(BVGetColumn(U,j,&u));
789 192 : PetscCall(VecPointwiseMult(u,u,svd->omega));
790 192 : PetscCall(BVRestoreColumn(U,j,&u));
791 192 : PetscCall(BVNormColumn(U,j,NORM_2,&nrm));
792 192 : PetscCall(BVScaleColumn(U,j,1.0/PetscAbs(nrm)));
793 192 : svd->sign[j] = PetscSign(nrm);
794 192 : PetscCall(BVNormColumn(V,j,NORM_2,&nrm));
795 192 : PetscCall(BVScaleColumn(V,j,1.0/nrm));
796 192 : j++;
797 : }
798 7 : PetscCall(VecDestroy(&x));
799 7 : PetscCall(VecDestroy(&x1));
800 7 : PetscCall(VecDestroy(&x2));
801 7 : PetscCall(VecDestroy(&xi));
802 7 : PetscCall(VecDestroy(&x1i));
803 7 : PetscCall(VecDestroy(&x2i));
804 7 : PetscFunctionReturn(PETSC_SUCCESS);
805 : }
806 :
807 38 : static PetscErrorCode SVDComputeVectors_Cyclic(SVD svd)
808 : {
809 38 : PetscFunctionBegin;
810 38 : switch (svd->problem_type) {
811 21 : case SVD_STANDARD:
812 21 : PetscCall(SVDComputeVectors_Cyclic_Standard(svd));
813 : break;
814 10 : case SVD_GENERALIZED:
815 10 : PetscCall(SVDComputeVectors_Cyclic_Generalized(svd));
816 : break;
817 7 : case SVD_HYPERBOLIC:
818 7 : PetscCall(SVDComputeVectors_Cyclic_Hyperbolic(svd));
819 : break;
820 0 : default:
821 0 : SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"Unknown singular value problem type");
822 : }
823 38 : PetscFunctionReturn(PETSC_SUCCESS);
824 : }
825 :
826 509 : static PetscErrorCode EPSMonitor_Cyclic(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
827 : {
828 509 : PetscInt i,j;
829 509 : SVD svd = (SVD)ctx;
830 509 : PetscScalar er,ei;
831 509 : PetscReal sigma;
832 509 : ST st;
833 :
834 509 : PetscFunctionBegin;
835 509 : nconv = 0;
836 509 : PetscCall(EPSGetST(eps,&st));
837 7759 : for (i=0,j=0;i<PetscMin(nest,svd->ncv);i++) {
838 7250 : er = eigr[i]; ei = eigi[i];
839 7250 : PetscCall(STBackTransform(st,1,&er,&ei));
840 7250 : PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL));
841 7250 : if (sigma>0.0) {
842 4713 : svd->sigma[j] = sigma;
843 4713 : svd->errest[j] = errest[i];
844 4713 : if (errest[i] && errest[i] < svd->tol) nconv++;
845 4713 : j++;
846 : }
847 : }
848 509 : nest = j;
849 509 : PetscCall(SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest));
850 509 : PetscFunctionReturn(PETSC_SUCCESS);
851 : }
852 :
853 30 : static PetscErrorCode SVDSetFromOptions_Cyclic(SVD svd,PetscOptionItems *PetscOptionsObject)
854 : {
855 30 : PetscBool set,val;
856 30 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
857 30 : ST st;
858 :
859 30 : PetscFunctionBegin;
860 30 : PetscOptionsHeadBegin(PetscOptionsObject,"SVD Cyclic Options");
861 :
862 30 : PetscCall(PetscOptionsBool("-svd_cyclic_explicitmatrix","Use cyclic explicit matrix","SVDCyclicSetExplicitMatrix",cyclic->explicitmatrix,&val,&set));
863 30 : if (set) PetscCall(SVDCyclicSetExplicitMatrix(svd,val));
864 :
865 30 : PetscOptionsHeadEnd();
866 :
867 30 : if (!cyclic->eps) PetscCall(SVDCyclicGetEPS(svd,&cyclic->eps));
868 30 : if (!cyclic->explicitmatrix && !cyclic->usereps) {
869 : /* use as default an ST with shell matrix and Jacobi */
870 18 : PetscCall(EPSGetST(cyclic->eps,&st));
871 18 : PetscCall(STSetMatMode(st,ST_MATMODE_SHELL));
872 : }
873 30 : PetscCall(EPSSetFromOptions(cyclic->eps));
874 30 : PetscFunctionReturn(PETSC_SUCCESS);
875 : }
876 :
877 19 : static PetscErrorCode SVDCyclicSetExplicitMatrix_Cyclic(SVD svd,PetscBool explicitmat)
878 : {
879 19 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
880 :
881 19 : PetscFunctionBegin;
882 19 : if (cyclic->explicitmatrix != explicitmat) {
883 12 : cyclic->explicitmatrix = explicitmat;
884 12 : svd->state = SVD_STATE_INITIAL;
885 : }
886 19 : PetscFunctionReturn(PETSC_SUCCESS);
887 : }
888 :
889 : /*@
890 : SVDCyclicSetExplicitMatrix - Indicate if the eigensolver operator
891 : H(A) = [ 0 A ; A^T 0 ] must be computed explicitly.
892 :
893 : Logically Collective
894 :
895 : Input Parameters:
896 : + svd - singular value solver
897 : - explicitmat - boolean flag indicating if H(A) is built explicitly
898 :
899 : Options Database Key:
900 : . -svd_cyclic_explicitmatrix <boolean> - Indicates the boolean flag
901 :
902 : Level: advanced
903 :
904 : .seealso: SVDCyclicGetExplicitMatrix()
905 : @*/
906 19 : PetscErrorCode SVDCyclicSetExplicitMatrix(SVD svd,PetscBool explicitmat)
907 : {
908 19 : PetscFunctionBegin;
909 19 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
910 57 : PetscValidLogicalCollectiveBool(svd,explicitmat,2);
911 19 : PetscTryMethod(svd,"SVDCyclicSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
912 19 : PetscFunctionReturn(PETSC_SUCCESS);
913 : }
914 :
915 1 : static PetscErrorCode SVDCyclicGetExplicitMatrix_Cyclic(SVD svd,PetscBool *explicitmat)
916 : {
917 1 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
918 :
919 1 : PetscFunctionBegin;
920 1 : *explicitmat = cyclic->explicitmatrix;
921 1 : PetscFunctionReturn(PETSC_SUCCESS);
922 : }
923 :
924 : /*@
925 : SVDCyclicGetExplicitMatrix - Returns the flag indicating if H(A) is built explicitly.
926 :
927 : Not Collective
928 :
929 : Input Parameter:
930 : . svd - singular value solver
931 :
932 : Output Parameter:
933 : . explicitmat - the mode flag
934 :
935 : Level: advanced
936 :
937 : .seealso: SVDCyclicSetExplicitMatrix()
938 : @*/
939 1 : PetscErrorCode SVDCyclicGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
940 : {
941 1 : PetscFunctionBegin;
942 1 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
943 1 : PetscAssertPointer(explicitmat,2);
944 1 : PetscUseMethod(svd,"SVDCyclicGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
945 1 : PetscFunctionReturn(PETSC_SUCCESS);
946 : }
947 :
948 1 : static PetscErrorCode SVDCyclicSetEPS_Cyclic(SVD svd,EPS eps)
949 : {
950 1 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
951 :
952 1 : PetscFunctionBegin;
953 1 : PetscCall(PetscObjectReference((PetscObject)eps));
954 1 : PetscCall(EPSDestroy(&cyclic->eps));
955 1 : cyclic->eps = eps;
956 1 : cyclic->usereps = PETSC_TRUE;
957 1 : svd->state = SVD_STATE_INITIAL;
958 1 : PetscFunctionReturn(PETSC_SUCCESS);
959 : }
960 :
961 : /*@
962 : SVDCyclicSetEPS - Associate an eigensolver object (EPS) to the
963 : singular value solver.
964 :
965 : Collective
966 :
967 : Input Parameters:
968 : + svd - singular value solver
969 : - eps - the eigensolver object
970 :
971 : Level: advanced
972 :
973 : .seealso: SVDCyclicGetEPS()
974 : @*/
975 1 : PetscErrorCode SVDCyclicSetEPS(SVD svd,EPS eps)
976 : {
977 1 : PetscFunctionBegin;
978 1 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
979 1 : PetscValidHeaderSpecific(eps,EPS_CLASSID,2);
980 1 : PetscCheckSameComm(svd,1,eps,2);
981 1 : PetscTryMethod(svd,"SVDCyclicSetEPS_C",(SVD,EPS),(svd,eps));
982 1 : PetscFunctionReturn(PETSC_SUCCESS);
983 : }
984 :
985 31 : static PetscErrorCode SVDCyclicGetEPS_Cyclic(SVD svd,EPS *eps)
986 : {
987 31 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
988 :
989 31 : PetscFunctionBegin;
990 31 : if (!cyclic->eps) {
991 31 : PetscCall(EPSCreate(PetscObjectComm((PetscObject)svd),&cyclic->eps));
992 31 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)cyclic->eps,(PetscObject)svd,1));
993 31 : PetscCall(EPSSetOptionsPrefix(cyclic->eps,((PetscObject)svd)->prefix));
994 31 : PetscCall(EPSAppendOptionsPrefix(cyclic->eps,"svd_cyclic_"));
995 31 : PetscCall(PetscObjectSetOptions((PetscObject)cyclic->eps,((PetscObject)svd)->options));
996 31 : PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL));
997 31 : PetscCall(EPSMonitorSet(cyclic->eps,EPSMonitor_Cyclic,svd,NULL));
998 : }
999 31 : *eps = cyclic->eps;
1000 31 : PetscFunctionReturn(PETSC_SUCCESS);
1001 : }
1002 :
1003 : /*@
1004 : SVDCyclicGetEPS - Retrieve the eigensolver object (EPS) associated
1005 : to the singular value solver.
1006 :
1007 : Collective
1008 :
1009 : Input Parameter:
1010 : . svd - singular value solver
1011 :
1012 : Output Parameter:
1013 : . eps - the eigensolver object
1014 :
1015 : Level: advanced
1016 :
1017 : .seealso: SVDCyclicSetEPS()
1018 : @*/
1019 31 : PetscErrorCode SVDCyclicGetEPS(SVD svd,EPS *eps)
1020 : {
1021 31 : PetscFunctionBegin;
1022 31 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
1023 31 : PetscAssertPointer(eps,2);
1024 31 : PetscUseMethod(svd,"SVDCyclicGetEPS_C",(SVD,EPS*),(svd,eps));
1025 31 : PetscFunctionReturn(PETSC_SUCCESS);
1026 : }
1027 :
1028 1 : static PetscErrorCode SVDView_Cyclic(SVD svd,PetscViewer viewer)
1029 : {
1030 1 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
1031 1 : PetscBool isascii;
1032 :
1033 1 : PetscFunctionBegin;
1034 1 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1035 1 : if (isascii) {
1036 1 : if (!cyclic->eps) PetscCall(SVDCyclicGetEPS(svd,&cyclic->eps));
1037 1 : PetscCall(PetscViewerASCIIPrintf(viewer," %s matrix\n",cyclic->explicitmatrix?"explicit":"implicit"));
1038 1 : PetscCall(PetscViewerASCIIPushTab(viewer));
1039 1 : PetscCall(EPSView(cyclic->eps,viewer));
1040 1 : PetscCall(PetscViewerASCIIPopTab(viewer));
1041 : }
1042 1 : PetscFunctionReturn(PETSC_SUCCESS);
1043 : }
1044 :
1045 33 : static PetscErrorCode SVDReset_Cyclic(SVD svd)
1046 : {
1047 33 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
1048 :
1049 33 : PetscFunctionBegin;
1050 33 : PetscCall(EPSReset(cyclic->eps));
1051 33 : PetscCall(MatDestroy(&cyclic->C));
1052 33 : PetscCall(MatDestroy(&cyclic->D));
1053 33 : PetscFunctionReturn(PETSC_SUCCESS);
1054 : }
1055 :
1056 32 : static PetscErrorCode SVDDestroy_Cyclic(SVD svd)
1057 : {
1058 32 : SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
1059 :
1060 32 : PetscFunctionBegin;
1061 32 : PetscCall(EPSDestroy(&cyclic->eps));
1062 32 : PetscCall(PetscFree(svd->data));
1063 32 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",NULL));
1064 32 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",NULL));
1065 32 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",NULL));
1066 32 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",NULL));
1067 32 : PetscFunctionReturn(PETSC_SUCCESS);
1068 : }
1069 :
1070 32 : SLEPC_EXTERN PetscErrorCode SVDCreate_Cyclic(SVD svd)
1071 : {
1072 32 : SVD_CYCLIC *cyclic;
1073 :
1074 32 : PetscFunctionBegin;
1075 32 : PetscCall(PetscNew(&cyclic));
1076 32 : svd->data = (void*)cyclic;
1077 32 : svd->ops->solve = SVDSolve_Cyclic;
1078 32 : svd->ops->solveg = SVDSolve_Cyclic;
1079 32 : svd->ops->solveh = SVDSolve_Cyclic;
1080 32 : svd->ops->setup = SVDSetUp_Cyclic;
1081 32 : svd->ops->setfromoptions = SVDSetFromOptions_Cyclic;
1082 32 : svd->ops->destroy = SVDDestroy_Cyclic;
1083 32 : svd->ops->reset = SVDReset_Cyclic;
1084 32 : svd->ops->view = SVDView_Cyclic;
1085 32 : svd->ops->computevectors = SVDComputeVectors_Cyclic;
1086 32 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",SVDCyclicSetEPS_Cyclic));
1087 32 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",SVDCyclicGetEPS_Cyclic));
1088 32 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",SVDCyclicSetExplicitMatrix_Cyclic));
1089 32 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",SVDCyclicGetExplicitMatrix_Cyclic));
1090 32 : PetscFunctionReturn(PETSC_SUCCESS);
1091 : }
|
