diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh index e02ed16a996f6795cac5a7ccbb4936772328bb0f..91c2bf158c738f231214363d45b976c57cd81711 100644 --- a/dune/microstructure/CorrectorComputer.hh +++ b/dune/microstructure/CorrectorComputer.hh @@ -817,16 +817,16 @@ public: unsigned int Solvertype = parameterSet_.get<unsigned int>("Solvertype", 1); unsigned int Solver_verbosity = parameterSet_.get<unsigned int>("Solver_verbosity", 2); - // printmatrix(std::cout, stiffnessMatrix_, "StiffnessMatrix", "--"); + x_1_.resize(load_alpha1_.size()); + x_2_.resize(load_alpha2_.size()); + x_3_.resize(load_alpha3_.size()); - // --- set initial values for solver - x_1_ = load_alpha1_; - x_2_ = load_alpha2_; - x_3_ = load_alpha3_; + // printmatrix(std::cout, stiffnessMatrix_, "StiffnessMatrix", "--"); if (parameterSet_.get<bool>("printMicroOutput", false)) std::cout << "start corrector solver..." << std::endl; Dune::Timer SolverTimer; + //////////////////////////////////////////////////////////////////////////////////// if (Solvertype==1) // CHOLMOD - SOLVER { if (parameterSet_.get<bool>("printMicroOutput", false)) @@ -837,14 +837,33 @@ public: std::cout << "WARNING: Using a direct solver with " << basis_.dimension() << " degrees of freedom may be not feasible or slow." << std::endl; Dune::Solvers::CholmodSolver<MatrixCT,VectorCT> solver; - solver.setProblem(stiffnessMatrix_,x_1_,load_alpha1_); - // solver.preprocess(); + + // Since we have the same stiffness-matrix for three different right-hand sides, + // we first use the factorization of CHOLMOD once and use it to repeatedly solve for multiple right-hand sides. + // In my measurements this reduced the runtime to about ~30% compared to using 'setProblem' multiple times. + + // set first solution vector + solver.setSolutionVector(x_1_); + // Factorize the matrix + solver.setMatrix(stiffnessMatrix_); + solver.factorize(); + // set first right-hand side + solver.setRhs(load_alpha1_); + // solve for first vector solver.solve(); - solver.setProblem(stiffnessMatrix_,x_2_,load_alpha2_); - // solver.preprocess(); + + // set second solution vector + solver.setSolutionVector(x_2_); + // set second right-hand side + solver.setRhs(load_alpha2_); + // solve for second vector solver.solve(); - solver.setProblem(stiffnessMatrix_,x_3_,load_alpha3_); - // solver.preprocess(); + + // set third solution vector + solver.setSolutionVector(x_3_); + // set third right-hand side + solver.setRhs(load_alpha3_); + // solve for third vector solver.solve(); // log_ << "Solver-type used: " <<" CHOLMOD-Solver" << std::endl; } @@ -875,6 +894,12 @@ public: { if (parameterSet_.get<bool>("printMicroOutput", false)) std::cout << "------------ GMRES - Solver ------------" << std::endl; + + // --- set initial values for solver (only for the iterativer solvers) + x_1_ = load_alpha1_; + x_2_ = load_alpha2_; + x_3_ = load_alpha3_; + // Turn the matrix into a linear operator Dune::MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_); @@ -911,6 +936,12 @@ public: { if (parameterSet_.get<bool>("printMicroOutput", false)) std::cout << "------------ CG - Solver ------------" << std::endl; + + // --- set initial values for solver (only for the iterativer solvers) + x_1_ = load_alpha1_; + x_2_ = load_alpha2_; + x_3_ = load_alpha3_; + Dune::MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_); // Sequential incomplete LU decomposition as the preconditioner @@ -950,6 +981,7 @@ public: if (parameterSet_.get<bool>("printMicroOutput", false)) std::cout << "------------ QR - Solver ------------" << std::endl; // log_ << "solveLinearSystems: We use QR solver.\n"; + //TODO install suitesparse Dune::SPQR<MatrixCT> sPQR(stiffnessMatrix_); // sPQR.setVerbosity(1);