Skip to content

Commit ada8440

Browse files
cmdougcmdoug
authored
Bug fix for getting left eigenvectors without breaking SVD (FreeFem#350)
1 parent ad4769d commit ada8440

1 file changed

Lines changed: 69 additions & 70 deletions

File tree

‎plugin/mpi/SLEPc-code.hpp‎

Lines changed: 69 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,7 @@ 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+
SVDSetInitialSpaces(svd, othern, otherbasis, n, basis);
389388
SVDSolve(svd);
390389
}
391390
else if(std::is_same<SType, NEP>::value) {
@@ -401,9 +400,9 @@ AnyType eigensolver<Type, K, SType>::E_eigensolver::operator()(Stack stack) cons
401400
for(int i = 0; i < n; ++i)
402401
VecDestroy(&basis[i]);
403402
delete [] basis;
404-
for(int i = 0; i < nl; ++i)
405-
VecDestroy(&lbasis[i]);
406-
delete [] lbasis;
403+
for(int i = 0; i < othern; ++i)
404+
VecDestroy(&otherbasis[i]);
405+
delete [] otherbasis;
407406
PetscInt nconv;
408407
if(std::is_same<SType, EPS>::value)
409408
EPSGetConverged(eps, &nconv);
@@ -413,79 +412,79 @@ AnyType eigensolver<Type, K, SType>::E_eigensolver::operator()(Stack stack) cons
413412
NEPGetConverged(nep, &nconv);
414413
else
415414
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])))) {
415+
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])))) {
417416
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;
418417
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;
418+
KNM<K>* eigenarray = nargs[4] ? GetAny<KNM<K>*>((*nargs[4])(stack)) : nullptr;
419+
KNM<K>* otherarray = nargs[8] ? GetAny<KNM<K>*>((*nargs[8])(stack)) : nullptr;
421420
if(eigenvalues)
422421
eigenvalues->resize(nconv);
423-
if(rvectors && !isType)
424-
rvectors->resize(nconv);
425-
if(lvectors && !isType)
426-
lvectors->resize(nconv);
422+
if(eigenvectors && !isType)
423+
eigenvectors->resize(nconv);
424+
if(othervectors && !isType)
425+
othervectors->resize(nconv);
427426
if(errorestimate)
428427
errorestimate->resize(nconv);
429-
if(rarray)
430-
rarray->resize(m, nconv);
428+
if(eigenarray)
429+
eigenarray->resize(m, nconv);
431430
Vec xr, xi, yr, yi;
432431
PetscInt n, nr;
433-
if(rvectors || rarray || lvectors || larray) {
432+
if(eigenvectors || eigenarray || othervectors || otherarray) {
434433
MatCreateVecs(ptA->_petsc, &xi, &xr);
435434
VecGetLocalSize(xi, &nr);
436-
if(larray)
437-
larray->resize(nr, nconv);
435+
if(otherarray)
436+
otherarray->resize(nr, nconv);
438437
VecGetLocalSize(xr, &n);
439438
} else xr = xi = NULL;
440-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)
439+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)
441440
MatCreateVecs(ptA->_petsc, &yi, &yr);
442441
else yr = yi = NULL;
443442
for(PetscInt i = 0; i < nconv; ++i) {
444443
PetscScalar kr, ki = 0;
445444
PetscReal sigma;
446445
PetscReal errest;
447446
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)
447+
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);
448+
if((othervectors || otherarray) && isTwoSided)
450449
EPSGetLeftEigenvector(eps, i, yr, yi);
451450
if(errorestimate)
452451
EPSGetErrorEstimate(eps, i, &errest);
453452
}
454453
else if(std::is_same<SType, SVD>::value)
455454
SVDGetSingularTriplet(svd, i, &sigma, xr, xi);
456455
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);
456+
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);
458457
if(errorestimate)
459458
NEPGetErrorEstimate(nep, i, &errest);
460459
}
461460
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);
461+
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);
463462
if(errorestimate)
464463
PEPGetErrorEstimate(pep, i, &errest);
465464
}
466-
if(rvectors || rarray || lvectors || larray) {
465+
if(eigenvectors || eigenarray || othervectors || otherarray) {
467466
PetscScalar* tmpr;
468467
PetscScalar* tmpi;
469468
PetscScalar* tmp2r;
470469
PetscScalar* tmp2i;
471470
VecGetArray(xr, &tmpr);
472-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)
471+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)
473472
VecGetArray(yr, &tmp2r);
474473
K* pt, *pti;
475474
if(!std::is_same<SType, SVD>::value) {
476475
if(std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value) {
477476
VecGetArray(xi, &tmpi);
478477
pt = new K[n];
479478
copy(pt, n, tmpr, tmpi);
480-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided) {
479+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided) {
481480
VecGetArray(yi, &tmp2i);
482-
pti = new K[nr];
483-
copy(pti, nr, tmp2r, tmp2i);
481+
pti = new K[n];
482+
copy(pti, n, tmp2r, tmp2i);
484483
}
485484
}
486485
else {
487486
pt = reinterpret_cast<K*>(tmpr);
488-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)
487+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)
489488
pti = reinterpret_cast<K*>(tmp2r);
490489
}
491490
}
@@ -502,13 +501,13 @@ AnyType eigensolver<Type, K, SType>::E_eigensolver::operator()(Stack stack) cons
502501
ptA->_A->HPDDM::template Subdomain<PetscScalar>::exchange(static_cast<K*>(cpy));
503502
else
504503
ptA->_exchange[0]->HPDDM::template Subdomain<PetscScalar>::exchange(static_cast<K*>(cpy));
505-
if(rvectors)
506-
rvectors->set(i, cpy);
507-
if(rarray && !codeA) {
504+
if(eigenvectors)
505+
eigenvectors->set(i, cpy);
506+
if(eigenarray && !codeA) {
508507
KN<K> cpy(m, pt);
509-
(*rarray)(':', i) = cpy;
508+
(*eigenarray)(':', i) = cpy;
510509
}
511-
if(std::is_same<SType, SVD>::value || (std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)) {
510+
if(std::is_same<SType, SVD>::value || (std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)) {
512511
KN<K> cpy(ptA->_cnum && ptA->_exchange[1] ? ptA->_exchange[1]->getDof() : (ptA->_A ? ptA->_A->getDof() : 0));
513512
cpy = K(0.0);
514513
if(ptA->_cnum && ptA->_exchange[1]) {
@@ -521,35 +520,35 @@ AnyType eigensolver<Type, K, SType>::E_eigensolver::operator()(Stack stack) cons
521520
}
522521
else
523522
ffassert(0);
524-
if(lvectors)
525-
lvectors->set(i, cpy);
526-
if(larray && !codeA) {
523+
if(othervectors)
524+
othervectors->set(i, cpy);
525+
if(otherarray && !codeA) {
527526
KN<K> cpy(nr, pti);
528-
(*larray)(':', i) = cpy;
527+
(*otherarray)(':', i) = cpy;
529528
}
530529
}
531530
}
532531
if(codeA || isType || !ptA->_A) {
533-
if(rarray) {
532+
if(eigenarray) {
534533
KN<K> cpy(m, pt);
535-
(*rarray)(':', i) = cpy;
534+
(*eigenarray)(':', i) = cpy;
536535
}
537-
if(larray) {
536+
if(otherarray) {
538537
KN<K> cpy(nr, pti);
539-
(*larray)(':', i) = cpy;
538+
(*otherarray)(':', i) = cpy;
540539
}
541540
}
542541
if(!std::is_same<SType, SVD>::value && std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value) {
543542
delete [] pt;
544-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)
543+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)
545544
delete [] pti;
546545
}
547546
if(std::is_same<SType, SVD>::value || (std::is_same<PetscScalar, double>::value && std::is_same<K, std::complex<double>>::value)) {
548547
VecRestoreArray(xi, &tmpi);
549-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)
548+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)
550549
VecRestoreArray(yi, &tmp2i);
551550
}
552-
if(std::is_same<SType, EPS>::value && (lvectors || larray) && isTwoSided)
551+
if(std::is_same<SType, EPS>::value && (othervectors || otherarray) && isTwoSided)
553552
VecRestoreArray(yr, &tmp2r);
554553
VecRestoreArray(xr, &tmpr);
555554
}
@@ -562,7 +561,7 @@ AnyType eigensolver<Type, K, SType>::E_eigensolver::operator()(Stack stack) cons
562561
if(errorestimate)
563562
errorestimate->operator[](i) = errest;
564563
}
565-
if(rvectors || rarray || lvectors || larray) {
564+
if(eigenvectors || eigenarray || othervectors || otherarray) {
566565
VecDestroy(&xr);
567566
VecDestroy(&xi);
568567
VecDestroy(&yr);

0 commit comments

Comments
 (0)