Skip to content

Commit 45060b0

Browse files
cmdougcmdoug
committed
Bug fix for getting left eigenvectors without breaking SVD
1 parent ad4769d commit 45060b0

1 file changed

Lines changed: 74 additions & 70 deletions

File tree

‎plugin/mpi/SLEPc-code.hpp‎

Lines changed: 74 additions & 70 deletions
Original file line numberDiff line numberDiff line change
@@ -147,12 +147,12 @@ basicAC_F0::name_and_type eigensolver<Type, K, SType>::E_eigensolver::name_param
147147
{"sparams", &typeid(std::string*)},
148148
{"prefix", &typeid(std::string*)},
149149
{"values", &typeid(KN<typename std::conditional<!std::is_same<SType, SVD>::value, K, PetscReal>::type>*)},
150-
{!std::is_same<SType, SVD>::value ? "vectors" : "rvectors", &typeid(FEbaseArrayKn<K>*)},
151-
{!std::is_same<SType, SVD>::value ? "array" : "rarray", &typeid(KNM<K>*)},
150+
{!std::is_same<SType, SVD>::value ? "vectors" : "lvectors", &typeid(FEbaseArrayKn<K>*)},
151+
{!std::is_same<SType, SVD>::value ? "array" : "larray", &typeid(KNM<K>*)},
152152
{"fields", &typeid(KN<double>*)},
153153
{"names", &typeid(KN<String>*)},
154-
{"lvectors", &typeid(FEbaseArrayKn<K>*)},
155-
{"larray", &typeid(KNM<K>*)},
154+
{!std::is_same<SType, SVD>::value ? "lvectors" : "rvectors", &typeid(FEbaseArrayKn<K>*)},
155+
{!std::is_same<SType, SVD>::value ? "larray" : "rarray", &typeid(KNM<K>*)},
156156
{"deflation", &typeid(KNM<PetscScalar>*)},
157157
{"errorestimate", &typeid(KN<PetscReal>*)},
158158
{"schurPreconditioner", &typeid(KN<Matrice_Creuse<upscaled_type<PetscScalar>>>*)},
@@ -324,50 +324,50 @@ AnyType eigensolver<Type, K, SType>::E_eigensolver::operator()(Stack stack) cons
324324
PEPSetFromOptions(pep);
325325
PEPSetUp(pep);
326326
}
327-
FEbaseArrayKn<K>* rvectors = nargs[3] ? GetAny<FEbaseArrayKn<K>*>((*nargs[3])(stack)) : nullptr;
327+
FEbaseArrayKn<K>* eigenvectors = nargs[3] ? GetAny<FEbaseArrayKn<K>*>((*nargs[3])(stack)) : nullptr;
328328
Vec* basis = nullptr;
329329
PetscInt n = 0;
330-
if(rvectors) {
330+
if(eigenvectors) {
331331
ffassert(!isType);
332-
if(rvectors->N > 0 && rvectors->get(0) && rvectors->get(0)->n > 0) {
333-
n = rvectors->N;
332+
if(eigenvectors->N > 0 && eigenvectors->get(0) && eigenvectors->get(0)->n > 0) {
333+
n = eigenvectors->N;
334334
basis = new Vec[n];
335335
for(int i = 0; i < n; ++i) {
336336
MatCreateVecs(ptA->_petsc, &basis[i], NULL);
337337
PetscScalar* pt;
338338
VecGetArray(basis[i], &pt);
339339
if(!(std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value))
340-
distributedVec(ptA->_num, ptA->_first, ptA->_last, static_cast<K*>(*(rvectors->get(i))), pt, rvectors->get(i)->n);
340+
distributedVec(ptA->_num, ptA->_first, ptA->_last, static_cast<K*>(*(eigenvectors->get(i))), pt, eigenvectors->get(i)->n);
341341
VecRestoreArray(basis[i], &pt);
342342
}
343343
}
344-
rvectors->resize(0);
344+
eigenvectors->resize(0);
345345
}
346-
FEbaseArrayKn<K>* lvectors = nargs[7] ? GetAny<FEbaseArrayKn<K>*>((*nargs[7])(stack)) : nullptr;
347-
Vec* lbasis = nullptr;
348-
PetscInt nl = 0;
349-
if(lvectors) {
346+
FEbaseArrayKn<K>* othervectors = nargs[7] ? GetAny<FEbaseArrayKn<K>*>((*nargs[7])(stack)) : nullptr;
347+
Vec* otherbasis = nullptr;
348+
PetscInt othern = 0;
349+
if(othervectors) {
350350
ffassert(!isType);
351-
if(lvectors->N > 0 && lvectors->get(0) && lvectors->get(0)->n > 0) {
352-
nl = lvectors->N;
353-
lbasis = new Vec[nl];
354-
for(int i = 0; i < nl; ++i) {
355-
MatCreateVecs(ptA->_petsc, &lbasis[i], NULL);
351+
if(othervectors->N > 0 && othervectors->get(0) && othervectors->get(0)->n > 0) {
352+
othern = othervectors->N;
353+
otherbasis = new Vec[othern];
354+
for(int i = 0; i < othern; ++i) {
355+
MatCreateVecs(ptA->_petsc, &otherbasis[i], NULL);
356356
PetscScalar* pt;
357-
VecGetArray(lbasis[i], &pt);
357+
VecGetArray(otherbasis[i], &pt);
358358
if(!(std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value))
359-
distributedVec(ptA->_num, ptA->_first, ptA->_last, static_cast<K*>(*(lvectors->get(i))), pt, lvectors->get(i)->n);
360-
VecRestoreArray(lbasis[i], &pt);
359+
distributedVec(ptA->_num, ptA->_first, ptA->_last, static_cast<K*>(*(othervectors->get(i))), pt, othervectors->get(i)->n);
360+
VecRestoreArray(otherbasis[i], &pt);
361361
}
362362
}
363-
lvectors->resize(0);
363+
othervectors->resize(0);
364364
}
365365
if(std::is_same<SType, EPS>::value) {
366366
EPSGetTwoSided(eps, &isTwoSided);
367367
if(n)
368368
EPSSetInitialSpace(eps, n, basis);
369-
if(nl && isTwoSided)
370-
EPSSetLeftInitialSpace(eps, nl, lbasis);
369+
if(othern && isTwoSided)
370+
EPSSetLeftInitialSpace(eps, othern, otherbasis);
371371
KNM<PetscScalar>* ptDeflation = nargs[9] ? GetAny<KNM<PetscScalar>*>((*nargs[9])(stack)) : NULL;
372372
if(ptDeflation && ptDeflation->M()) {
373373
PetscInt m;
@@ -384,8 +384,12 @@ AnyType eigensolver<Type, K, SType>::E_eigensolver::operator()(Stack stack) cons
384384
EPSSolve(eps);
385385
}
386386
else if(std::is_same<SType, SVD>::value) {
387-
if(n)
388-
SVDSetInitialSpaces(svd, n, basis, 0, NULL);
387+
if(n && !othern)
388+
SVDSetInitialSpaces(svd, 0, NULL, n, basis);
389+
else if(!n && othern)
390+
SVDSetInitialSpaces(svd, othern, otherbasis, 0, NULL);
391+
else if(n && othern)
392+
SVDSetInitialSpaces(svd, othern, otherbasis, n, basis);
389393
SVDSolve(svd);
390394
}
391395
else if(std::is_same<SType, NEP>::value) {
@@ -401,9 +405,9 @@ AnyType eigensolver<Type, K, SType>::E_eigensolver::operator()(Stack stack) cons
401405
for(int i = 0; i < n; ++i)
402406
VecDestroy(&basis[i]);
403407
delete [] basis;
404-
for(int i = 0; i < nl; ++i)
405-
VecDestroy(&lbasis[i]);
406-
delete [] lbasis;
408+
for(int i = 0; i < othern; ++i)
409+
VecDestroy(&otherbasis[i]);
410+
delete [] otherbasis;
407411
PetscInt nconv;
408412
if(std::is_same<SType, EPS>::value)
409413
EPSGetConverged(eps, &nconv);
@@ -413,79 +417,79 @@ AnyType eigensolver<Type, K, SType>::E_eigensolver::operator()(Stack stack) cons
413417
NEPGetConverged(nep, &nconv);
414418
else
415419
PEPGetConverged(pep, &nconv);
416-
if(nconv > 0 && ((nargs[2] || nargs[3] || nargs[4]) || ((std::is_same<SType, SVD>::value || std::is_same<SType, EPS>::value) && (nargs[7] || nargs[8])))) {
420+
if(nconv > 0 && ((nargs[2] || nargs[3] || nargs[4]) || ((std::is_same<SType, SVD>::value || (std::is_same<SType, EPS>::value && isTwoSided)) && (nargs[7] || nargs[8])))) {
417421
KN<typename std::conditional<!std::is_same<SType, SVD>::value, K, PetscReal>::type>* eigenvalues = nargs[2] ? GetAny<KN<typename std::conditional<!std::is_same<SType, SVD>::value, K, PetscReal>::type>*>((*nargs[2])(stack)) : nullptr;
418422
KN<PetscReal>* errorestimate = nargs[10] ? GetAny<KN<PetscReal>*>((*nargs[10])(stack)) : nullptr;
419-
KNM<K>* rarray = nargs[4] ? GetAny<KNM<K>*>((*nargs[4])(stack)) : nullptr;
420-
KNM<K>* larray = nargs[8] ? GetAny<KNM<K>*>((*nargs[8])(stack)) : nullptr;
423+
KNM<K>* eigenarray = nargs[4] ? GetAny<KNM<K>*>((*nargs[4])(stack)) : nullptr;
424+
KNM<K>* otherarray = nargs[8] ? GetAny<KNM<K>*>((*nargs[8])(stack)) : nullptr;
421425
if(eigenvalues)
422426
eigenvalues->resize(nconv);
423-
if(rvectors && !isType)
424-
rvectors->resize(nconv);
425-
if(lvectors && !isType)
426-
lvectors->resize(nconv);
427+
if(eigenvectors && !isType)
428+
eigenvectors->resize(nconv);
429+
if(othervectors && !isType)
430+
othervectors->resize(nconv);
427431
if(errorestimate)
428432
errorestimate->resize(nconv);
429-
if(rarray)
430-
rarray->resize(m, nconv);
433+
if(eigenarray)
434+
eigenarray->resize(m, nconv);
431435
Vec xr, xi, yr, yi;
432436
PetscInt n, nr;
433-
if(rvectors || rarray || lvectors || larray) {
437+
if(eigenvectors || eigenarray || othervectors || otherarray) {
434438
MatCreateVecs(ptA->_petsc, &xi, &xr);
435439
VecGetLocalSize(xi, &nr);
436-
if(larray)
437-
larray->resize(nr, nconv);
440+
if(otherarray)
441+
otherarray->resize(nr, nconv);
438442
VecGetLocalSize(xr, &n);
439443
} else xr = xi = NULL;
440-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)
444+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)
441445
MatCreateVecs(ptA->_petsc, &yi, &yr);
442446
else yr = yi = NULL;
443447
for(PetscInt i = 0; i < nconv; ++i) {
444448
PetscScalar kr, ki = 0;
445449
PetscReal sigma;
446450
PetscReal errest;
447451
if(std::is_same<SType, EPS>::value) {
448-
EPSGetEigenpair(eps, i, &kr, &ki, (rvectors || rarray) ? xr : NULL, (rvectors || rarray) && std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value ? xi : NULL);
449-
if((lvectors || larray) && isTwoSided)
452+
EPSGetEigenpair(eps, i, &kr, &ki, (eigenvectors || eigenarray) ? xr : NULL, (eigenvectors || eigenarray) && std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value ? xi : NULL);
453+
if((othervectors || otherarray) && isTwoSided)
450454
EPSGetLeftEigenvector(eps, i, yr, yi);
451455
if(errorestimate)
452456
EPSGetErrorEstimate(eps, i, &errest);
453457
}
454458
else if(std::is_same<SType, SVD>::value)
455459
SVDGetSingularTriplet(svd, i, &sigma, xr, xi);
456460
else if(std::is_same<SType, NEP>::value) {
457-
NEPGetEigenpair(nep, i, &kr, &ki, (rvectors || rarray) ? xr : NULL, (rvectors || rarray) && std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value ? xi : NULL);
461+
NEPGetEigenpair(nep, i, &kr, &ki, (eigenvectors || eigenarray) ? xr : NULL, (eigenvectors || eigenarray) && std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value ? xi : NULL);
458462
if(errorestimate)
459463
NEPGetErrorEstimate(nep, i, &errest);
460464
}
461465
else {
462-
PEPGetEigenpair(pep, i, &kr, &ki, (rvectors || rarray) ? xr : NULL, (rvectors || rarray) && std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value ? xi : NULL);
466+
PEPGetEigenpair(pep, i, &kr, &ki, (eigenvectors || eigenarray) ? xr : NULL, (eigenvectors || eigenarray) && std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value ? xi : NULL);
463467
if(errorestimate)
464468
PEPGetErrorEstimate(pep, i, &errest);
465469
}
466-
if(rvectors || rarray || lvectors || larray) {
470+
if(eigenvectors || eigenarray || othervectors || otherarray) {
467471
PetscScalar* tmpr;
468472
PetscScalar* tmpi;
469473
PetscScalar* tmp2r;
470474
PetscScalar* tmp2i;
471475
VecGetArray(xr, &tmpr);
472-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)
476+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)
473477
VecGetArray(yr, &tmp2r);
474478
K* pt, *pti;
475479
if(!std::is_same<SType, SVD>::value) {
476480
if(std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value) {
477481
VecGetArray(xi, &tmpi);
478482
pt = new K[n];
479483
copy(pt, n, tmpr, tmpi);
480-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided) {
484+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided) {
481485
VecGetArray(yi, &tmp2i);
482-
pti = new K[nr];
483-
copy(pti, nr, tmp2r, tmp2i);
486+
pti = new K[n];
487+
copy(pti, n, tmp2r, tmp2i);
484488
}
485489
}
486490
else {
487491
pt = reinterpret_cast<K*>(tmpr);
488-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)
492+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)
489493
pti = reinterpret_cast<K*>(tmp2r);
490494
}
491495
}
@@ -502,13 +506,13 @@ AnyType eigensolver<Type, K, SType>::E_eigensolver::operator()(Stack stack) cons
502506
ptA->_A->HPDDM::template Subdomain<PetscScalar>::exchange(static_cast<K*>(cpy));
503507
else
504508
ptA->_exchange[0]->HPDDM::template Subdomain<PetscScalar>::exchange(static_cast<K*>(cpy));
505-
if(rvectors)
506-
rvectors->set(i, cpy);
507-
if(rarray && !codeA) {
509+
if(eigenvectors)
510+
eigenvectors->set(i, cpy);
511+
if(eigenarray && !codeA) {
508512
KN<K> cpy(m, pt);
509-
(*rarray)(':', i) = cpy;
513+
(*eigenarray)(':', i) = cpy;
510514
}
511-
if(std::is_same<SType, SVD>::value || (std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)) {
515+
if(std::is_same<SType, SVD>::value || (std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)) {
512516
KN<K> cpy(ptA->_cnum && ptA->_exchange[1] ? ptA->_exchange[1]->getDof() : (ptA->_A ? ptA->_A->getDof() : 0));
513517
cpy = K(0.0);
514518
if(ptA->_cnum && ptA->_exchange[1]) {
@@ -521,35 +525,35 @@ AnyType eigensolver<Type, K, SType>::E_eigensolver::operator()(Stack stack) cons
521525
}
522526
else
523527
ffassert(0);
524-
if(lvectors)
525-
lvectors->set(i, cpy);
526-
if(larray && !codeA) {
528+
if(othervectors)
529+
othervectors->set(i, cpy);
530+
if(otherarray && !codeA) {
527531
KN<K> cpy(nr, pti);
528-
(*larray)(':', i) = cpy;
532+
(*otherarray)(':', i) = cpy;
529533
}
530534
}
531535
}
532536
if(codeA || isType || !ptA->_A) {
533-
if(rarray) {
537+
if(eigenarray) {
534538
KN<K> cpy(m, pt);
535-
(*rarray)(':', i) = cpy;
539+
(*eigenarray)(':', i) = cpy;
536540
}
537-
if(larray) {
541+
if(otherarray) {
538542
KN<K> cpy(nr, pti);
539-
(*larray)(':', i) = cpy;
543+
(*otherarray)(':', i) = cpy;
540544
}
541545
}
542546
if(!std::is_same<SType, SVD>::value && std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value) {
543547
delete [] pt;
544-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)
548+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)
545549
delete [] pti;
546550
}
547551
if(std::is_same<SType, SVD>::value || (std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value)) {
548552
VecRestoreArray(xi, &tmpi);
549-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)
553+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)
550554
VecRestoreArray(yi, &tmp2i);
551555
}
552-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)
556+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)
553557
VecRestoreArray(yr, &tmp2r);
554558
VecRestoreArray(xr, &tmpr);
555559
}
@@ -562,7 +566,7 @@ AnyType eigensolver<Type, K, SType>::E_eigensolver::operator()(Stack stack) cons
562566
if(errorestimate)
563567
errorestimate->operator[](i) = errest;
564568
}
565-
if(rvectors || rarray || lvectors || larray) {
569+
if(eigenvectors || eigenarray || othervectors || otherarray) {
566570
VecDestroy(&xr);
567571
VecDestroy(&xi);
568572
VecDestroy(&yr);

0 commit comments

Comments
 (0)