ug4
Loading...
Searching...
No Matches
pinvit.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 * Author: Martin Rupp
4 *
5 * This file is part of UG4.
6 *
7 * UG4 is free software: you can redistribute it and/or modify it under the
8 * terms of the GNU Lesser General Public License version 3 (as published by the
9 * Free Software Foundation) with the following additional attribution
10 * requirements (according to LGPL/GPL v3 §7):
11 *
12 * (1) The following notice must be displayed in the Appropriate Legal Notices
13 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 *
15 * (2) The following notice must be displayed at a prominent place in the
16 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 *
18 * (3) The following bibliography is recommended for citation and must be
19 * preserved in all covered files:
20 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 * parallel geometric multigrid solver on hierarchically distributed grids.
22 * Computing and visualization in science 16, 4 (2013), 151-164"
23 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 * flexible software system for simulating pde based models on high performance
25 * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 *
27 * This program is distributed in the hope that it will be useful,
28 * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 * GNU Lesser General Public License for more details.
31 */
32
33#ifndef __H__UG__LIB_ALGEBRA__PINVIT_H__
34#define __H__UG__LIB_ALGEBRA__PINVIT_H__
35
36#include <complex>
37#include <vector>
38#include <string>
41#include "additional_math.h"
42
44
49#include "common/progress.h"
50
52
53#include "smart_ptr_vector.h"
54
55// constructors
56namespace ug{
57
58
101template<typename TAlgebra>
103 : public DebugWritingObject<TAlgebra>
104{
105public:
106// Algebra type
107 typedef TAlgebra algebra_type;
108
109// Vector type
110 typedef typename TAlgebra::vector_type vector_type;
111 typedef typename TAlgebra::matrix_type matrix_type;
112 typedef typename vector_type::value_type vector_value_type;
113
116
117
118private:
122
123 //ILinearOperator<vector_type,vector_type>* m_pA;
124 //ILinearOperator<vector_type,vector_type>* m_pB;
126
129
131
134
138 size_t m_iPINVIT;
139 std::vector<double> lambda;
140 std::vector<bool> vbDirichlet;
141
143
145
146 std::vector<double> freq_change;
151
162
163public:
164 std::string tostring()
165 {
166 return config_string();
167 }
168 virtual std::string config_string() const
169 {
170 std::stringstream ss;
171 ss << "PINVIT Eigensolver by Martin Rupp / G-CSC 2013-2015." <<
172 "\n MaxIterations = " << m_maxIterations <<
173 "\n Precision = " << m_dPrecision << (m_bRelativePrecision ? " (relative) " : " (absolute) ") <<
174 "\n MinimumDefectToCalcCorrection = " << m_dMinimumDefectToCalcCorrection <<
175 "\n Number of EV = " << px.size() <<
176 "\n PINVIT = " << m_iPINVIT << " = ";
177 if(m_iPINVIT == 1)
178 ss << "Preconditioned Inverse Block Iteration [Neymeyr]";
179 else if(m_iPINVIT==2) ss << "Preconditioned Block Gradient Method";
180 else if(m_iPINVIT==3) ss << "LOBPCG (locally optimal block preconditioned gradient) [Knyazew]";
181 else if(m_iPINVIT>=4) ss << "gerneralized method PINVIT-method looking back " << m_iPINVIT-1 << " ev approximations.";
182 ss << "\n Preconditioner: ";
183 if(m_spPrecond.valid()) ss << ConfigShift(m_spPrecond->config_string()) << "\n";
184 else ss << " NOT SET!\n";
185
187 {
188 ss << "\n\tAdditionaly storing " << m_additionalEigenvectorsToKeep << " eigenvectors";
189 }
191 ss << "\n\tUsing additional Corrections.";
192 if(m_bLaplacian)
193 ss << "\n\tUsing LAPLACIAN!";
194
195 ss << "\n";
196 return ss.str();
197
198 }
199
229
230 void set_abort_on_frequencies_converged_linear_elasticity(double freq_precision, double density)
231 {
235 m_dFreqPrecision = freq_precision;
236 }
237
239 {
241 }
247 {
248 px.push_back(vec);
249 }
250
257 {
258 m_spPrecond = precond;
259 }
260
262 {
263 m_pA = loA;
264 // get dirichlet nodes
266 m_pA.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
267 matrix_type &A = pOp->get_matrix();
268 vbDirichlet.resize(A.num_rows());
269
270 if(m_bLaplacian)
271 {
272 bool bFirst=false;
273 for(size_t i=0; i<A.num_rows(); i++)
274 {
275 //vbDirichlet[i] = A.is_isolated(i);
276 if(bFirst && A.is_isolated(i))
277 {
278 bFirst = false;
279 continue;
280 }
281
282 for(typename matrix_type::row_iterator it = A.begin_row(i); it != A.end_row(i); ++it)
283 {
284 it.value() = -1.0;
285 }
286 A(i, i) = A.num_connections(i)*1.0;
287
288 }
289
290 for(size_t i=0; i<A.num_rows(); i++)
291 vbDirichlet[i] = false;
292 }
293 else
294 {
295 for(size_t i=0; i<A.num_rows(); i++)
296 vbDirichlet[i] = A.is_isolated(i);
297 }
298
299 }
300
302 {
303 m_pB = &B;
304 }
305
306 void set_max_iterations(size_t maxIterations)
307 {
308 m_maxIterations = maxIterations;
309 }
310
311 void set_precision(double precision)
312 {
313 m_dPrecision = precision;
315 m_bRelativePrecision = false;
316 }
317
322
327
332
340 void set_pinvit(size_t iPINVIT)
341 {
342 m_iPINVIT = iPINVIT;
343 UG_ASSERT(iPINVIT >=1 && iPINVIT <= 3, "i has to be >= 1 and <= 3, but is " << iPINVIT);
344 }
345
346 size_t num_eigenvalues()
347 {
348 return px.size();
349 }
350
351 double get_eigenvalue(size_t i)
352 {
353 if(lambda.size() != px.size()) return 0.0;
354 return lambda[i];
355 }
356
357 SmartPtr<vector_type> get_eigenvector(size_t i)
358 {
359 return px[i];
360 }
361
362 void set_print_eigenvalues_and_defect(bool b)
363 {
364 m_bPrintEigenvaluesAndDefect = b;
365 }
366
367 void set_print_projected_eigenvectors(bool b)
368 {
369 m_bPrintProjectedEigenvectors = b;
370 }
371 void set_print_projected_eigenvalues(bool b)
372 {
373 m_bPrintProjectedEigenvalues = b;
374 }
375
376 void set_print_projected_eigenproblem(bool b)
377 {
378 m_bPrintProjectedEigenproblem = b;
379 }
380
381 void set_print_used_testvectors(bool b)
382 {
383 m_bPrintUsedTestvectors = b;
384 }
385
386 void print_projected_eigenvectors(DenseMatrix<VariableArray2<double> > &r_ev, std::vector<std::complex<double> > &r_lambda, std::vector<std::string> &vTestVectorDescription)
387 {
388 UG_LOG("evs: \n");
389 for(size_t c=0; c < r_ev.num_cols(); c++)
390 {
391 UG_LOG("eigenvalue [" << c << "] = " << r_lambda[c] << ", vector:\n");
392 std::vector<double> tmpEV(r_ev.num_rows());
393 for(size_t r=0; r<r_ev.num_rows(); r++)
394 tmpEV[r] = r_ev(r, c);
395 std::vector<size_t> s = GetSortedIndices(tmpEV, absCompare);
396
397
398 for(size_t i=0; i<r_ev.num_rows(); i++)
399 {
400 //if(r_ev.num_rows() > 3 && abs(r_ev(s[i], c))/abs(r_ev(s[r_ev.num_rows()-1], c)) < 0.01) continue;
401 size_t j = s[i];
402 UG_LOG(std::setw(14) << r_ev(j, c) << " " << vTestVectorDescription[j] << "\n");
403 }
404 UG_LOG("\n\n");
405 }
406
407
408 }
409
410 SmartPtr<vector_type> create_approximation_vector()
411 {
412 SmartPtr<vector_type> t = px(0).clone_without_values();
413 t->set(0.0);
414 return t;
415 }
416
417 void set_store_defects(bool b)
418 {
419 m_bStoreDefects = b;
420 }
421
422 void set_store_lambdas(bool b)
423 {
424 m_bStoreLambdas = b;
425 }
426
427 bool m_bStoreDefects;
428 bool m_bStoreLambdas;
429 size_t m_iteration;
430 std::vector<std::vector<double> > m_defects;
431 std::vector<std::vector<double> > m_lambdas;
432
433 size_t get_iterations()
434 {
435 return m_iteration;
436 }
437
438 double get_defect(size_t nev, size_t it)
439 {
440 UG_COND_THROW(m_bStoreDefects == false, "Not storing defects. Use set_store_defects(true).");
441 UG_COND_THROW(nev > m_defects.size(), nev << " > #EV" << m_defects.size());
442 UG_COND_THROW(it > m_defects[nev].size(), "Iteration " << it << " not available.");
443 return m_defects[nev][it];
444 }
445
446 double get_lambda(size_t nev, size_t it)
447 {
448 UG_COND_THROW(m_bStoreLambdas == false, "Not storing lambdas. Use set_store_lambdas(true).");
449 UG_COND_THROW(nev > m_lambdas.size(), nev << " > #EV" << m_lambdas.size());
450 UG_COND_THROW(it > m_lambdas[nev].size(), "Iteration " << it << " not available.");
451
452 return m_lambdas[nev][it];
453 }
454
459 int apply()
460 {
461 PINVIT_PROFILE_FUNC();
462 UG_LOG(config_string() << "\nInitializing... ");
463
464
465 size_t nEigenvalues = px.size();
466 DenseMatrix<VariableArray2<double> > rA;
467 DenseMatrix<VariableArray2<double> > rB;
468 DenseMatrix<VariableArray2<double> > r_ev;
469 std::vector<std::complex<double> > r_lambda;
470 std::vector<std::string> vCorrectionName;
471
472 SmartPtrVector<vector_type> pTestVectors;
473 SmartPtrVector<vector_type> vAdditional;
474 SmartPtrVector<vector_type> vCorr;
475 SmartPtrVector<vector_type> vOldX;
476
477 std::vector<double> vDefectNorm(nEigenvalues, m_dPrecision*10);
478 std::vector<double> oldXnorm(nEigenvalues);
479 std::vector<bool> bConverged(nEigenvalues, false);
480
481 std::vector<std::string> vTestVectorDescription;
482
483 SmartPtr<vector_type> spDefect;
484
485 try{
486 spDefect = create_approximation_vector();
487
488 m_currentAdditionalEigenvectors = 0;
489
490 //size_t size = px[0]->size();
491 /*
492 ParallelMatrix<SparseMatrix<double> > B;
493 B.resize(size, size);
494 for(size_t i=0; i<size; i++)
495 B(i,i) = 0.00390625;
496 B.set_storage_type(PST_ADDITIVE);*/
497
498 for(size_t i=0; i<m_additionalEigenvectorsToKeep; i++)
499 vAdditional.push_back(create_approximation_vector());
500
501 lambda.resize(nEigenvalues);
502 freq_change.resize(nEigenvalues);
503 for(size_t i=0; i<nEigenvalues; i++)
504 {
505 vCorr.push_back(create_approximation_vector() );
506 if(m_iPINVIT >= 3)
507 vOldX.push_back(create_approximation_vector() );
508 }
509
510 vCorrectionName.resize(nEigenvalues);
511 }UG_CATCH_THROW("PINVIT::apply init vectors failed.");
512
513 UG_LOG("done.\n");
514 UG_LOG("PINVIT: Initializing preconditioner... ");
515
516 try{
517 m_spPrecond->init(m_pA);
518 }UG_CATCH_THROW("PINVIT::apply init preconditioner.");
519
520 UG_LOG("done.\n");
521
522#ifdef UG_PARALLEL
523 pcl::ProcessCommunicator pc;
524#endif
525
526 m_iteration=0;
527 m_defects.clear();
528 if(m_bStoreDefects)
529 m_defects.resize(nEigenvalues);
530 m_lambdas.clear();
531 if(m_bStoreLambdas)
532 m_lambdas.resize(nEigenvalues);
533
534 for(m_iteration=0; m_iteration<m_maxIterations; m_iteration++)
535 {
536 try{
537
538 // UG_LOG("normalize...\n");
539 // 0. normalize
540 normalize_approximations();
541
542 // 2. compute rayleigh quotient, residuum, apply preconditioner, compute corrections norm
543 size_t nrofconverged=0;
544
545
546 size_t currentAdditionalCorrections = 0;
547 size_t numCorrections=0;
548
549 PROGRESS_START(prog, nEigenvalues, "PINVIT: compute corrections");
550 for(size_t i=0; i<nEigenvalues; i++)
551 {
552 PROGRESS_UPDATE(prog, i);
553 // UG_LOG("compute rayleigh " << i << "...\n");
554
555 double freq_old = 0;
556
557 if(m_bPrintFrequencies_linear_elasticity){
558 freq_old = 1.0/(2.0*PI)*sqrt(lambda[i]/m_dDensity_linear_elasticity);
559 }
560
561 compute_rayleigh_and_defect(px(i), lambda[i], *spDefect, vDefectNorm[i]);
562
563 if(m_bPrintFrequencies_linear_elasticity){
564 freq_change[i] = freq_old-1.0/(2.0*PI)*sqrt(lambda[i]/m_dDensity_linear_elasticity);
565 }
566
567 if((m_iteration != 0 && m_bRelativePrecision && vDefectNorm[i]/lambda[i] < m_dPrecision)
568 || (!m_bRelativePrecision && vDefectNorm[i] < m_dPrecision))
569 bConverged[i] = true;
570 else
571 bConverged[i] = false;
572
573 if(bConverged[i])
574 {
575 nrofconverged++;
576 if(!m_bRelativePrecision &&
577 m_bUseAdditionalCorrections && currentAdditionalCorrections < m_currentAdditionalEigenvectors)
578 {
579 double additionalLambda, additionalDefect;
580 compute_rayleigh_and_defect(vAdditional(currentAdditionalCorrections),
581 additionalLambda, *spDefect, additionalDefect);
582 currentAdditionalCorrections++;
583 if(additionalDefect > m_dPrecision && additionalDefect < 1e5)
584 {
585 calculate_correction(*vCorr[numCorrections++], *spDefect);
586 vCorrectionName[numCorrections-1] = std::string("additional correction ") + ToString(currentAdditionalCorrections-1);
587 }
588 }
589 }
590 else
591 {
592 // UG_LOG("Calculating correction from defect " << i << "\n");
593 vCorrectionName[numCorrections] = std::string("correction ") + ToString(i);
594 calculate_correction(*vCorr[numCorrections++], *spDefect);
595 }
596
597 if(m_bStoreDefects)
598 m_defects[i].push_back(vDefectNorm[i]);
599 if(m_bStoreLambdas)
600 m_lambdas[i].push_back(lambda[i]);
601 }
602 PROGRESS_FINISH(prog);
603
604
605
606
607 // output
608 if(m_bPrintEigenvaluesAndDefect)
609 print_eigenvalues_and_defect(m_iteration, vDefectNorm, oldXnorm, lambda, bConverged);
610
611 if(nrofconverged==nEigenvalues)
612 {
613 UG_LOG("all eigenvectors converged\n");
614 return true;
615 }
616
617 if(m_bAbortOnFreqConverged_linear_elasticity){
618 bool s = true;
619 for(unsigned k = 0; k < freq_change.size(); ++k){
620 if(abs(freq_change[k]) >= m_dFreqPrecision){
621 s = false;
622 break;
623 }
624 }
625
626 if(s){
627 UG_LOG("all eigenvectors converged according to frequency change\n");
628 return true;
629 }
630 }
631
632 // 5. add Testvectors
633 //UG_LOG("5. add Testvectors\nEigenvalues");
634
635 get_testvectors(m_iteration, vCorr, vCorrectionName, numCorrections, vOldX, vAdditional, pTestVectors, vTestVectorDescription, bConverged);
636
637 for(size_t i=0; i<nEigenvalues; i++)
638 {
639 write_debug(m_iteration, i, px(i), *spDefect, vCorr(i), bConverged[i]);
640 if(vOldX.size() > i) write_debug_old(m_iteration, i, vOldX(i));
641 }
642
643 /*for(size_t i=0; i<vTestVectorDescription.size(); i++)
644 { UG_LOG(vTestVectorDescription[i] << "\nEigenvalues"); } */
645
646 // 5. compute reduced Matrices rA, rB
647
648 if(pTestVectors.size() < px.size()+1)
649 {
650 UG_LOG("ERROR: Projected Spaces has Dimension " << pTestVectors.size() << ", but we need at least " << px.size()+1 << " to make progress (for " << px.size() << " eigenvalues)\n");
651 return false;
652 }
653
654 get_projected_eigenvalue_problem(rA, rB, pTestVectors, vTestVectorDescription);
655
656 // 6. solve reduced eigenvalue problem
657 size_t iNrOfTestVectors = pTestVectors.size();
658 r_ev.resize(iNrOfTestVectors, iNrOfTestVectors);
659 r_lambda.resize(iNrOfTestVectors);
660
661 // solve rA x = lambda rB, --> r_ev, r_lambda
662 GeneralizedEigenvalueProblemComplex(rA, r_ev, r_lambda, rB, true);
663
664 if(m_bPrintProjectedEigenvalues)
665 print_projected_eigenvectors(r_ev, r_lambda, vTestVectorDescription);
666
667
668 if(m_bDebugCalcProjectedEigenvalues)
669 debug_calc_projected_eigenvalues(r_ev, pTestVectors, m_iteration, false);
670 set_new_approximations_and_save_old(r_ev, pTestVectors, vOldX, vAdditional);
671
672 assert_real_positive(r_lambda);
673 }UG_CATCH_THROW("PINVIT::apply iteration " << m_iteration << " failed.");
674 }
675
676 UG_LOG("not converged after" << m_maxIterations << " steps.\n");
677 return false;
678 }
679
680 // only use this type of vector creation if you don't use the preconditioner
681 void init_nonprecond_vector(vector_type &t)
682 {
683 CloneVector(t, px(0));
684 t.set(0.0);
685 }
693 void debug_calc_projected_eigenvalues(DenseMatrix<VariableArray2<double> > &r_ev,
694 SmartPtrVector<vector_type> &pTestVectors, int iteration, bool bWrite)
695 {
696 PROFILE_FUNC_GROUP("debug");
697 try{
698 if(debug_writer() == SPNULL) return;
699
700 vector_type t, defect;
701 init_nonprecond_vector(t);
702 init_nonprecond_vector(defect);
703
704 #ifdef UG_PARALLEL
705 for(size_t c=0; c<pTestVectors.size(); c++)
706 pTestVectors[c]->change_storage_type(PST_CONSISTENT);
707 #endif
708
709 for(size_t r=0; r<pTestVectors.size(); r++)
710 {
711 #ifdef UG_PARALLEL
712 t.set_storage_type(PST_CONSISTENT);
713 defect.set_storage_type(PST_ADDITIVE);
714 #endif
715 for(size_t i=0; i<t.size(); i++)
716 {
717 t[i] = 0.0;
718 if(!vbDirichlet[i])
719 {
720 for(size_t c=0; c<pTestVectors.size(); c++)
721 t[i] += r_ev(c, r) * (*pTestVectors[c])[i];
722 }
723 }
724
725
726 t *= 1.0/B_norm(t);
727
728 #ifdef UG_PARALLEL
729 t.set_storage_type(PST_CONSISTENT);
730 defect.set_storage_type(PST_ADDITIVE);
731 #endif
732 m_pA->apply(defect, t);
733
734 double lambda = t.dotprod(defect);
735
736 if(m_pB)
737 {
738 #ifdef UG_PARALLEL
739 defect.change_storage_type(PST_ADDITIVE);
740 t.change_storage_type(PST_CONSISTENT);
741 #endif
742 MatMultAdd(defect, 1.0, defect, -lambda, *m_pB, t);
743 }
744 else
745 VecScaleAdd(defect, 1.0, defect, -lambda, t);
746
747
748
749 UG_LOG(r << " lambda: " << std::setw(14) << lambda << " defect: " << std::setw(14) << defect.norm() << "\n");
750
751 if(bWrite)
752 write_debug(t, ("pinvit_it_" + ToString(iteration) + "_pev_" + ToString(r)).c_str());
753 }
754 }UG_CATCH_THROW("PINVIT::debug_calc_projected_eigenvalues failed.");
755 }
756
757
758 void set_new_approximations_and_save_old(DenseMatrix<VariableArray2<double> > &r_ev, SmartPtrVector<vector_type> &pTestVectors, SmartPtrVector<vector_type> &vOldX,
759 SmartPtrVector<vector_type> &vAdditional)
760 {
761 PINVIT_PROFILE_FUNC();
762 try{
763 #ifdef UG_PARALLEL
764 for(size_t i=0; i<pTestVectors.size(); i++)
765 pTestVectors[i]->change_storage_type(PST_CONSISTENT);
766 for(size_t i=0; i<px.size(); i++)
767 px[i]->change_storage_type(PST_CONSISTENT);
768
769 for(size_t i=0; i<vOldX.size(); i++)
770 vOldX[i]->set_storage_type(PST_CONSISTENT);
771
772 for(size_t i=0; i<vAdditional.size(); i++)
773 vAdditional[i]->set_storage_type(PST_CONSISTENT);
774 #endif
775 // assume r_lambda is sorted
776
777 size_t size = px[0]->size();
778 THROW_IF_NOT_EQUAL(pTestVectors.size(), r_ev.num_rows());
779
780 size_t numCalc = std::min(px.size()+m_additionalEigenvectorsToKeep, r_ev.num_cols()); //px.size()
781 m_currentAdditionalEigenvectors = numCalc-px.size();
782 // UG_LOG("current aditional eigenvectors: " << m_currentAdditionalEigenvectors << "\n");
783
784 std::vector<vector_value_type> x_tmp(numCalc);
785 for(size_t i=0; i<size; i++)
786 {
787 if(vbDirichlet[i]) continue;
788
789 // since vx can be part of the Testvectors, temporary safe result in x_tmp.
790 for(size_t r=0; r < numCalc; r++)
791 {
792 x_tmp[r] = 0.0;
793 for(size_t c=0; c<pTestVectors.size(); c++)
794 x_tmp[r] += r_ev(c, r) * (*pTestVectors[c])[i];
795 }
796
797 // now overwrite
798 for(size_t r=0; r<px.size(); r++)
799 {
800 // save old
801 /*for(int k=0; k<m_nrOfOld; k++)
802 vOldX[k][r][i] = vOldX[k][r][i];
803 vOldX[0][r][i] = (px[r])[i];*/
804 if(vOldX.size()>0)
805 vOldX(r) [i] = px(r) [i];
806 // store new
807 px(r) [i] = x_tmp[r];
808 }
809
810 // store additional
811 for(size_t r=px.size(); r<numCalc; r++)
812 vAdditional(r-px.size()) [i] = x_tmp[r];
813
814 }
815 }UG_CATCH_THROW("PINVIT::set_new_approximations_and_save_old failed.");
816 }
817
818 void assert_real_positive(const std::vector<std::complex<double> > &r_lambda)
819 {
820 PINVIT_PROFILE_FUNC();
821 try{
822 //#define UG_ASSERT2(cond, txt) {if(!cond) { UG_LOG("Warning: " << txt << "\n"); }}
823
824 for(size_t i=0; i<r_lambda.size(); i++) // px.size() or r_lambda.size()
825 {
826 if(r_lambda[i].imag() >= 1e-8)
827 {UG_LOG("WARNING: eigenvalue " << i << " is imaginary (" << r_lambda[i] << ")\n"); }
828 if(r_lambda[i].real() <= 1e-8)
829 {UG_LOG("WARNING: eigenvalues " << i << "<= 0\n");}
830 }
831
832 if(m_bPrintProjectedEigenvalues && m_bPrintProjectedEigenvectors == false)
833 {
834 for(size_t i=0; i<r_lambda.size(); i++)
835 UG_LOG(i << ".: " << r_lambda[i] << "\n");
836 UG_LOG("\n");
837 }
838
839 for(size_t i=0; i<r_lambda.size(); i++) // px.size() or r_lambda.size()
840 {
841 UG_ASSERT(r_lambda[i].imag() < 1e-4 && r_lambda[i].imag() > -1e-4, "eigenvalue " << i << " = " << r_lambda[i] << " is imaginary (" << r_lambda[i] << ")\n");
842 UG_ASSERT(r_lambda[i].real() > -1e-4, "eigenvalues " << i << " = " << r_lambda[i] << " < 0\n");
843 }
844 }UG_CATCH_THROW("PINVIT::assert_real_positive failed.");
845
846 }
847
848 void get_max_deflection_of_a_mode(vector_type& maxDeflect, vector_type& statSol,
849 const vector_type& eigenVec, const matrix_type& massMat)
850 {
851 PINVIT_PROFILE_FUNC();
852#ifdef UG_PARALLEL
853 statSol.set_storage_type(PST_CONSISTENT);
854 maxDeflect.set_storage_type(PST_CONSISTENT);
855#endif
856
857 SmartPtr<vector_type> Bv = statSol.clone_without_values();
858 SmartPtr<vector_type> vwBv = statSol.clone_without_values();
859
860 // compute Bv = massMat * eigenVec;
861 massMat.axpy(*Bv, 0.0, eigenVec, 1.0, eigenVec);
862
863 // vwBv = eigenVec * < statSol, Bv >
864 number scalarProd = 0.0;
865 for(size_t i = 0; i < statSol.size(); i++)
866 scalarProd += VecProd((*Bv)[i], statSol[i]);
867
868 for(size_t i = 0; i < statSol.size(); i++)
869 {
870 (*vwBv)[i] = eigenVec[i] * (-1.0) * scalarProd;
871 // maxDeflect += vwBv
872 maxDeflect[i] = maxDeflect[i] + (*vwBv)[i];
873 }
874
875 }
876
877private:
878 void write_debug(int iteration, int i, vector_type &x, vector_type &defect, vector_type &corr, bool bConverged)
879 {
880 PROFILE_FUNC_GROUP("debug");
881 write_debug(x, ("pinvit_it_" + ToString(iteration) + "_ev_" + ToString(i)).c_str());
882 write_debug(defect, ("pinvit_it_" + ToString(iteration) + "_defect_" + ToString(i)).c_str());
883 if(!bConverged)
884 write_debug(corr, ("pinvit_it_" + ToString(iteration) + "_corr_" + ToString(i)).c_str());
885 }
886
887 void write_debug_old(int iteration, int i, vector_type &oldX)
888 {
889 write_debug(oldX, ("pinvit_it_" + ToString(iteration) + "_old_" + ToString(i)).c_str());
890 }
891
892 double B_norm(vector_type &x)
893 {
894 PINVIT_PROFILE_FUNC();
895 if(m_pB != NULL)
896 return EnergyNorm(x, *m_pB);
897 else
898 return x.norm();
899 }
900
901 void normalize_approximations()
902 {
903 PINVIT_PROFILE_FUNC();
904 for(size_t i=0; i< px.size(); i++)
905 px(i) *= 1/ B_norm(px(i));
906 }
907
908 void calculate_correction(vector_type &corr, vector_type &defect)
909 {
910 PINVIT_PROFILE_FUNC();
911 try{
912 #ifdef UG_PARALLEL
913 defect.change_storage_type(PST_ADDITIVE);
914 corr.set_storage_type(PST_CONSISTENT);
915 #endif
916 corr *= 0.0;
917 // d. apply preconditioner
918 if(1)
919 {
920 m_spPrecond->apply(corr, defect);
921 }
922 else
923 {
924 m_spPrecond->apply_update_defect(corr, defect);
925 }
926 corr *= 1/ B_norm(corr);
927 #ifdef UG_PARALLEL
928 corr.change_storage_type(PST_UNIQUE);
929 #endif
930 }UG_CATCH_THROW("PINVIT::calculate_correction failed.");
931 }
941 void compute_rayleigh_and_defect(vector_type &x, double &lambda, vector_type &defect, double &defectNorm)
942 {
943 PINVIT_PROFILE_FUNC();
944 try{
945 // a. compute rayleigh quotients
946
947 #ifdef UG_PARALLEL
948 x.change_storage_type(PST_CONSISTENT);
949 defect.set_storage_type(PST_ADDITIVE);
950 #endif
951 m_pA->apply(defect, x);
952
953
954 #ifdef UG_PARALLEL
955 defect.change_storage_type(PST_UNIQUE);
956 x.change_storage_type(PST_UNIQUE);
957 #endif
958 lambda = x.dotprod(defect); // / <px[i], Bpx[i]> = 1.0.
959 //UG_LOG("lambda[" << i << "] = " << lambda << "\n");
960
961 // b. calculate residuum
962 // defect = A px[i] - lambda[i] B px[i]
963 if(m_pB)
964 {
965 #ifdef UG_PARALLEL
966 defect.change_storage_type(PST_ADDITIVE);
967 x.change_storage_type(PST_CONSISTENT);
968 #endif
969 MatMultAdd(defect, 1.0, defect, -lambda, *m_pB, x);
970
971 }
972 else
973 VecScaleAdd(defect, 1.0, defect, -lambda, x);
974
975 // c. check if converged
976
977
978 #ifdef UG_PARALLEL
979 defect.change_storage_type(PST_UNIQUE);
980 #endif
981 defectNorm = defect.norm();
982 }UG_CATCH_THROW("PINVIT::compute_rayleigh_and_defect failed.");
983 }
984
992 void print_eigenvalues_and_defect(int iteration, const std::vector<double> &vDefectNorm,
993 std::vector<double> &vOldDefectNorm, const std::vector<double> &vLambda, std::vector<bool> &bConverged)
994 {
995 PINVIT_PROFILE_FUNC();
996 UG_LOG("=====================================================================================\n");
997 UG_LOG("iteration " << iteration << "\n");
998
999 for(size_t i=0; i<vLambda.size(); i++)
1000 {
1001 UG_LOG(i << " lambda: " << std::setw(20) << vLambda[i]);
1002 if(m_bPrintFrequencies_linear_elasticity){
1003 UG_LOG(" freq: " << std::setw(20) << 1.0/(2.0*PI)*sqrt(vLambda[i]/m_dDensity_linear_elasticity) );
1004 UG_LOG(" change: " << std::setw(20) << freq_change[i]);
1005 }
1006
1007 UG_LOG(" defect: " << std::setw(20) << vDefectNorm[i]);
1008 if(iteration != 0) { UG_LOG(" reduction: " << std::setw(14) << vDefectNorm[i]/vOldDefectNorm[i]); }
1009 if(m_bPrintFrequencies_linear_elasticity && abs(freq_change[i]) < m_dFreqPrecision) { UG_LOG(" (freq. converged) "); }
1010 if(bConverged[i]) { UG_LOG(" (converged)"); }
1011 UG_LOG("\n");
1012 vOldDefectNorm[i] = vDefectNorm[i];
1013 }
1014 UG_LOG("\n");
1015 }
1016
1034 void get_testvectors(int iteration, SmartPtrVector<vector_type> &vCorr,
1035 std::vector<std::string> &vCorrectionName,
1036 size_t numCorrections, SmartPtrVector<vector_type> &vOldX, SmartPtrVector<vector_type> &vAdditional,
1037 SmartPtrVector<vector_type> &pTestVectors, std::vector<std::string> &vTestVectorDescription,
1038 std::vector<bool> &bConverged)
1039 {
1040 try{
1041 PINVIT_PROFILE_FUNC();
1042 pTestVectors.clear();
1043 vTestVectorDescription.clear();
1044 if(m_iPINVIT == 1)
1045 {
1046 UG_ASSERT(0, "todo");
1047 for(size_t i=0; i < px.size(); i++)
1048 {
1049 if(!bConverged[i])
1050 {
1051 VecScaleAdd(vCorr(i), -1.0, vCorr(i), 1.0, px(i));
1052 vTestVectorDescription.push_back(std::string("ev - vCorr [") + ToString(i) + std::string("]") );
1053 }
1054 }
1055 }
1056 else
1057 {
1058 for(size_t i=0; i < px.size(); i++)
1059 {
1060 if(IsFiniteAndNotTooBig(px(i)))
1061 {
1062 pTestVectors.push_back(px[i]);
1063 vTestVectorDescription.push_back(std::string("eigenvector [") + ToString(i) + std::string("]") );
1064
1065 if(!bConverged[i] && m_iPINVIT >= 3)
1066 {
1067 if(iteration != 0)
1068 {
1069 UG_ASSERT(vOldX.size() > i, "?");
1070 pTestVectors.push_back(vOldX[i]);
1071 vTestVectorDescription.push_back(std::string("old eigenvalue[") + ToString(i) + std::string("]") );
1072 }
1073
1074 /*if(iteration == 0)
1075 {
1076 for(size_t j=0; j<px[i]->size(); j++)
1077 vOldX[i][j] = (px[i])[j] * urand(-1.0, 1.0);
1078 }*/
1079 }
1080 }
1081 else
1082 { UG_LOG("WARNING ev [" << i << "] not normal or too big\n"); }
1083 }
1084
1085 for(size_t i=0; i<numCorrections; i++)
1086 {
1087 if(IsFiniteAndNotTooBig(vCorr(i)))
1088 {
1089 pTestVectors.push_back(vCorr[i]);
1090 vTestVectorDescription.push_back(vCorrectionName[i]);
1091 }
1092 else
1093 {
1094 UG_LOG("WARNING correction [" << i << "] not normal or too big\n");
1095 }
1096 }
1097
1098 }
1099 for(size_t i=0; i<m_currentAdditionalEigenvectors; i++)
1100 {
1101 if(IsFiniteAndNotTooBig(vAdditional(i)))
1102 {
1103 pTestVectors.push_back(vAdditional[i]);
1104 vTestVectorDescription.push_back(std::string("additional [") + ToString(i) + std::string("]") );
1105 }
1106 else
1107 { UG_LOG("WARNING additional [" << i << "] not normal or too big\n"); }
1108 }
1109 }UG_CATCH_THROW("PINVIT::get_testvectors failed.");
1110 }
1111
1112
1118 void get_linear_independent_rows(DenseMatrix<VariableArray2<double> > mat, std::vector<bool> &bLinearIndependent)
1119 {
1120 PINVIT_PROFILE_FUNC();
1121 // Remove linear depended vectors
1122 bLinearIndependent.resize(mat.num_rows(), true);
1123 for(size_t i=0; i<mat.num_rows(); i++)
1124 {
1125 for(size_t j=0; j<i; j++)
1126 {
1127 if(!bLinearIndependent[i]) continue;
1128 double val = mat(i, j)/mat(j,j);
1129 mat(i,j) = 0;
1130 for(size_t k=j+1; k<mat.num_rows(); k++)
1131 mat(i,k) -= val*mat(j, k);
1132 }
1133 if(mat(i,i) < m_linearDependentEps) bLinearIndependent[i] = false;
1134 else bLinearIndependent[i] = true;
1135 }
1136 }
1137
1143 template<typename T>
1144 void remove_unused(T &v, const std::vector<bool> vbUse)
1145 {
1146 PINVIT_PROFILE_FUNC();
1147 T tmp;
1148 tmp = v;
1149 v.clear();
1150 for(size_t i=0; i<tmp.size(); i++)
1151 if(vbUse[i])
1152 v.push_back(tmp[i]);
1153
1154 }
1155
1156 void print_used_testvectors(std::vector<std::string> &vTestVectorDescription, std::vector<bool> bUse)
1157 {
1158 UG_LOG("used testvectors:\n");
1159 for(size_t i=0; i<vTestVectorDescription.size(); i++)
1160 if(bUse[i]) { UG_LOG(vTestVectorDescription[i] << "\n"); }
1161 UG_LOG("unused testvectors:\n");
1162 for(size_t i=0; i<vTestVectorDescription.size(); i++)
1163 if(!bUse[i]) { UG_LOG(vTestVectorDescription[i] << "\n"); }
1164 }
1165
1176 void get_projected_eigenvalue_problem(DenseMatrix<VariableArray2<double> > &rA,
1177 DenseMatrix<VariableArray2<double> > &rB, SmartPtrVector<vector_type> &pTestVectors,
1178 std::vector<std::string> &vTestVectorDescription)
1179 {
1180 try{
1181 PINVIT_PROFILE_FUNC();
1182 // 1. calculate W as a subset of the testvectors so that those are linear B-independent
1183
1184 size_t iNrOfTestVectors = pTestVectors.size();
1185 rA.resize(iNrOfTestVectors, iNrOfTestVectors);
1186 rB.resize(iNrOfTestVectors, iNrOfTestVectors);
1187
1188 if(m_pB)
1189 MultiEnergyProd(*m_pB, pTestVectors, rB, iNrOfTestVectors);
1190 else
1191 MultiScalProd(pTestVectors, rB, iNrOfTestVectors);
1192
1193
1194 // Remove linear depended vectors
1195 std::vector<bool> bUse;
1196 get_linear_independent_rows(rB, bUse);
1197
1198
1199 if(m_bPrintUsedTestvectors)
1200 print_used_testvectors(vTestVectorDescription, bUse);
1201
1202 // save used testvectors
1203 remove_unused(pTestVectors, bUse);
1204 remove_unused(vTestVectorDescription, bUse);
1205
1206 iNrOfTestVectors = pTestVectors.size();
1207
1208 // 2. & 3. compute reduced Matrices rA, rB
1209 rA.resize(iNrOfTestVectors, iNrOfTestVectors);
1210 rB.resize(iNrOfTestVectors, iNrOfTestVectors);
1211
1212 // UG_LOG("iNrOfTestVectors = " << iNrOfTestVectors << "\n");
1213
1214 if(m_pB)
1215 MultiEnergyProd(*m_pB, pTestVectors, rB, iNrOfTestVectors);
1216 else
1217 MultiScalProd(pTestVectors, rB, iNrOfTestVectors);
1218
1219 MultiEnergyProd(*m_pA, pTestVectors, rA, iNrOfTestVectors);
1220
1221 if(m_bPrintProjectedEigenproblem)
1222 {
1223 PrintMaple(rA, "rA");
1224 PrintMaple(rB, "rB");
1225 }
1226 }UG_CATCH_THROW("PINVIT::get_projected_eigenvalue_problem failed.");
1227
1228 }
1229public:
1230 void set_laplacian(bool b)
1231 {
1232 m_bLaplacian = b;
1233 }
1234
1235 void set_relative_precision(double precision)
1236 {
1237 m_dPrecision = precision;
1238 m_dMinimumDefectToCalcCorrection = precision;
1239 m_bRelativePrecision = true;
1240 }
1241
1242};
1243
1244} // namespace ug
1245
1246
1247#endif // __H__UG__LIB_ALGEBRA__PINVIT_H__
Definition smart_pointer.h:108
Definition debug_writer.h:354
virtual void set_debug(SmartPtr< IDebugWriter< algebra_type > > spDebugWriter)
set debug writer
Definition debug_writer.h:384
void write_debug(const matrix_type &mat, const char *filename)
write debug output for a matrix (if debug writer set)
Definition debug_writer.h:399
SmartPtr< IDebugWriter< algebra_type > > debug_writer()
returns the debug writer
Definition debug_writer.h:391
describes a linear iterator
Definition linear_iterator.h:81
describes a linear mapping X->Y
Definition linear_operator.h:80
Definition pinvit.h:104
size_t m_currentAdditionalCorrections
Definition pinvit.h:161
void set_store_defects(bool b)
Definition pinvit.h:417
bool m_bRelativePrecision
Definition pinvit.h:137
bool m_bPrintFrequencies_linear_elasticity
Definition pinvit.h:148
SmartPtr< ILinearIterator< vector_type > > m_spPrecond
Definition pinvit.h:133
void set_debug_calc_projected_eigenvalues(bool b)
Definition pinvit.h:323
void set_new_approximations_and_save_old(DenseMatrix< VariableArray2< double > > &r_ev, SmartPtrVector< vector_type > &pTestVectors, SmartPtrVector< vector_type > &vOldX, SmartPtrVector< vector_type > &vAdditional)
Definition pinvit.h:758
void set_linear_operator_B(matrix_operator_type &B)
Definition pinvit.h:301
void set_use_additional_corrections(bool b)
Definition pinvit.h:328
size_t m_currentAdditionalEigenvectors
Definition pinvit.h:160
double m_dFreqPrecision
Definition pinvit.h:150
std::vector< double > freq_change
Definition pinvit.h:146
SmartPtr< ILinearOperator< vector_type > > m_pA
Definition pinvit.h:125
std::string tostring()
Definition pinvit.h:164
void set_abort_on_frequencies_converged_linear_elasticity(double freq_precision, double density)
Definition pinvit.h:230
void set_linear_dependent_eps(double eps)
Definition pinvit.h:238
matrix_operator_type * m_pB
Definition pinvit.h:128
TAlgebra::vector_type vector_type
Definition pinvit.h:110
void calculate_correction(vector_type &corr, vector_type &defect)
Definition pinvit.h:908
bool m_bPrintEigenvaluesAndDefect
Definition pinvit.h:152
bool m_bStoreDefects
Definition pinvit.h:427
int apply()
Definition pinvit.h:459
void add_vector(SmartPtr< vector_type > vec)
Definition pinvit.h:246
std::vector< bool > vbDirichlet
Definition pinvit.h:140
double m_dDensity_linear_elasticity
Definition pinvit.h:147
double m_linearDependentEps
Definition pinvit.h:142
void set_store_lambdas(bool b)
Definition pinvit.h:422
DebugWritingObject< TAlgebra > base_type
Base type.
Definition pinvit.h:115
double m_dPrecision
Definition pinvit.h:136
size_t m_iPINVIT
Definition pinvit.h:138
bool m_bUseAdditionalCorrections
Definition pinvit.h:157
size_t m_additionalEigenvectorsToKeep
Definition pinvit.h:159
IPreconditioner< TAlgebra >::matrix_operator_type matrix_operator_type
Definition pinvit.h:127
std::vector< double > lambda
Definition pinvit.h:139
bool m_bLaplacian
Definition pinvit.h:144
TAlgebra::matrix_type matrix_type
Definition pinvit.h:111
void set_preconditioner(SmartPtr< ILinearIterator< vector_type > > precond)
Definition pinvit.h:256
bool m_bStoreLambdas
Definition pinvit.h:428
size_t m_maxIterations
Definition pinvit.h:135
bool m_bPrintUsedTestvectors
Definition pinvit.h:156
vector_type::value_type vector_value_type
Definition pinvit.h:112
virtual std::string config_string() const
Definition pinvit.h:168
bool m_bPrintProjectedEigenproblem
Definition pinvit.h:155
void set_precision(double precision)
Definition pinvit.h:311
void debug_calc_projected_eigenvalues(DenseMatrix< VariableArray2< double > > &r_ev, SmartPtrVector< vector_type > &pTestVectors, int iteration, bool bWrite)
Definition pinvit.h:693
void set_pinvit(size_t iPINVIT)
Definition pinvit.h:340
void assert_real_positive(const std::vector< std::complex< double > > &r_lambda)
Definition pinvit.h:818
bool m_bPrintProjectedEigenvalues
Definition pinvit.h:153
void set_max_iterations(size_t maxIterations)
Definition pinvit.h:306
void get_projected_eigenvalue_problem(DenseMatrix< VariableArray2< double > > &rA, DenseMatrix< VariableArray2< double > > &rB, SmartPtrVector< vector_type > &pTestVectors, std::vector< std::string > &vTestVectorDescription)
Definition pinvit.h:1176
void set_additional_eigenvectors_to_keep(size_t i)
Definition pinvit.h:318
void set_linear_operator_A(SmartPtr< ILinearOperator< vector_type > > loA)
Definition pinvit.h:261
SmartPtrVector< vector_type > px
Definition pinvit.h:132
bool m_bPrintProjectedEigenvectors
Definition pinvit.h:154
bool m_bAbortOnFreqConverged_linear_elasticity
Definition pinvit.h:149
TAlgebra algebra_type
Definition pinvit.h:107
PINVIT()
Definition pinvit.h:200
void compute_rayleigh_and_defect(vector_type &x, double &lambda, vector_type &defect, double &defectNorm)
Definition pinvit.h:941
bool m_bDebugCalcProjectedEigenvalues
Definition pinvit.h:158
double m_dMinimumDefectToCalcCorrection
Definition pinvit.h:130
void get_testvectors(int iteration, SmartPtrVector< vector_type > &vCorr, std::vector< std::string > &vCorrectionName, size_t numCorrections, SmartPtrVector< vector_type > &vOldX, SmartPtrVector< vector_type > &vAdditional, SmartPtrVector< vector_type > &pTestVectors, std::vector< std::string > &vTestVectorDescription, std::vector< bool > &bConverged)
Definition pinvit.h:1034
Definition smart_ptr_vector.h:41
function util rates kinetic compute(ConvRateSetup)
virtual void init()
const NullSmartPtr SPNULL
The equivalent to NULL for smart pointers.
Definition smart_pointer.h:90
#define UG_ASSERT(expr, msg)
Definition assert.h:70
the ug namespace
string ConfigShift(string s)
Definition string_util.cpp:457
for(yylen=0;yystr[yylen];yylen++) continue