HermesCommon 1.0
mumps_solver.cpp
Go to the documentation of this file.
```00001 // This file is part of HermesCommon
00002 //
00003 // Copyright (c) 2009 hp-FEM group at the University of Nevada, Reno (UNR).
00005 //
00006 // Hermes2D is free software; you can redistribute it and/or modify
00007 // it under the terms of the GNU General Public License as published
00008 // by the Free Software Foundation; either version 2 of the License,
00009 // or (at your option) any later version.
00010 //
00011 // Hermes2D is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014 // GNU General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU General Public License
00017 // along with Hermes2D; if not, write to the Free Software
00018 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
00022 #include "config.h"
00023 #ifdef WITH_MUMPS
00024 #include "mumps_solver.h"
00025 #include "callstack.h"
00026
00027 namespace Hermes
00028 {
00029   namespace Algebra
00030   {
00031     extern "C"
00032     {
00033       extern void dmumps_c(DMUMPS_STRUC_C *mumps_param_ptr);
00034       extern void zmumps_c(ZMUMPS_STRUC_C *mumps_param_ptr);
00035     }
00036
00037 #define USE_COMM_WORLD  -987654
00038
00045     static int find_position(int *Ai, int Alen, int idx)
00046     {
00047       assert(idx >= 0);
00048
00049       register int lo = 0, hi = Alen - 1, mid;
00050
00051       while (1)
00052       {
00053         mid = (lo + hi) >> 1;
00054
00055         if(idx < Ai[mid]) hi = mid - 1;
00056         else if(idx > Ai[mid]) lo = mid + 1;
00057         else break;
00058
00060         // value to this position, return 0 when obtaining value there).
00061         if(lo > hi)
00062         {
00063           mid = -1;
00064           break;
00065         }
00066       }
00067       return mid;
00068     }
00069
00070     template<typename Scalar>
00071     MumpsMatrix<Scalar>::MumpsMatrix()
00072     {
00073       nnz = 0;
00074       this->size = 0;
00075       irn = NULL;
00076       jcn = NULL;
00077       Ax = NULL;
00078       Ap = NULL;
00079       Ai = NULL;
00080     }
00081
00082     template<typename Scalar>
00083     MumpsMatrix<Scalar>::~MumpsMatrix()
00084     {
00085       free();
00086     }
00087
00088     template<typename Scalar>
00089     void MumpsMatrix<Scalar>::alloc()
00090     {
00091       assert(this->pages != NULL);
00092
00093       // initialize the arrays Ap and Ai
00094       Ap = new unsigned int[this->size + 1];
00095       int aisize = this->get_num_indices();
00096       Ai = new int[aisize];
00097
00098       // sort the indices and remove duplicities, insert into Ai
00099       unsigned int i, pos = 0;
00100       for (i = 0; i < this->size; i++)
00101       {
00102         Ap[i] = pos;
00103         pos += sort_and_store_indices(this->pages[i], Ai + pos, Ai + aisize);
00104       }
00105       Ap[i] = pos;
00106
00107       delete [] this->pages;
00108       this->pages = NULL;
00109
00110       nnz = Ap[this->size];
00111
00112       Ax = new typename mumps_type<Scalar>::mumps_Scalar[nnz];
00113       memset(Ax, 0, sizeof(Scalar) * nnz);
00114
00115       irn = new int[nnz];
00116       jcn = new int[nnz];
00117       for (unsigned int i = 0; i < nnz; i++)
00118       {
00119         irn[i] = 1;
00120         jcn[i] = 1;
00121       }
00122     }
00123
00124     template<typename Scalar>
00125     void MumpsMatrix<Scalar>::free()
00126     {
00127       nnz = 0;
00128       delete [] Ap; Ap = NULL;
00129       delete [] Ai; Ai = NULL;
00130       delete [] Ax; Ax = NULL;
00131       delete [] irn; irn = NULL;
00132       delete [] jcn; jcn = NULL;
00133     }
00134
00135     inline double mumps_to_Scalar(double x)
00136     {
00137       return x;
00138     }
00139
00140     inline std::complex<double> mumps_to_Scalar(ZMUMPS_COMPLEX x)
00141     {
00142       return std::complex<double>(x.r, x.i);
00143     }
00144
00145     template<typename Scalar>
00146     Scalar MumpsMatrix<Scalar>::get(unsigned int m, unsigned int n)
00147     {
00148       // Find m-th row in the n-th column.
00149       int mid = find_position(Ai + Ap[n], Ap[n + 1] - Ap[n], m);
00150       // Return 0 if the entry has not been found.
00151       if(mid < 0) return 0.0;
00152       // Otherwise, add offset to the n-th column and return the value.
00153       if(mid >= 0) mid += Ap[n];
00154       return mumps_to_Scalar(Ax[mid]);
00155     }
00156
00157     template<typename Scalar>
00158     void MumpsMatrix<Scalar>::zero()
00159     {
00160       memset(Ax, 0, sizeof(Scalar) * Ap[this->size]);
00161     }
00162
00163     inline ZMUMPS_COMPLEX& operator +=(ZMUMPS_COMPLEX &a, std::complex<double> b)
00164     {
00165       a.r +=b.real();
00166       a.i +=b.imag();
00167       return a;
00168     }
00169
00170     template<typename Scalar>
00171     void MumpsMatrix<Scalar>::add(unsigned int m, unsigned int n, Scalar v)
00172     {
00173       //          produced an error in neutronics-2-group-adapt (although tutorial-07
00174       //          ran well).
00175       // Find m-th row in the n-th column.
00176       int pos = find_position(Ai + Ap[n], Ap[n + 1] - Ap[n], m);
00177       // Make sure we are adding to an existing non-zero entry.
00178       if(pos < 0)
00180       // Add offset to the n-th column.
00181       pos += Ap[n];
00183       Ax[pos] += v;
00184       irn[pos] = m + 1;  // MUMPS is indexing from 1
00185       jcn[pos] = n + 1;
00186     }
00187
00188     template<typename Scalar>
00189     void MumpsMatrix<Scalar>::add(unsigned int m, unsigned int n, Scalar **mat, int *rows, int *cols)
00190     {
00191       for (unsigned int i = 0; i < m; i++)       // rows
00192         for (unsigned int j = 0; j < n; j++)     // cols
00193           if(rows[i] >= 0 && cols[j] >= 0) // not Dir. dofs.
00195     }
00196
00198
00199     template<typename Scalar>
00201     {
00202       for (unsigned int i = 0; i < this->size; i++)
00203       {
00205       }
00206     };
00207
00210     template<typename Scalar>
00211     bool MumpsMatrix<Scalar>::dump(FILE *file, const char *var_name, EMatrixDumpFormat fmt, char* number_format)
00212     {
00213       // TODO
00214       switch (fmt)
00215       {
00216       case DF_PLAIN_ASCII:
00217         fprintf(file, "%d\n", this->size);
00218         fprintf(file, "%d\n", nnz);
00219         for (unsigned int i = 0; i < nnz; i++)
00220         {
00221           fprintf(file, "%d %d ", irn[i], jcn[i]);
00222           Hermes::Helpers::fprint_num(file, mumps_to_Scalar(Ax[i]));
00223           fprintf(file, "\n");
00224         }
00225         return true;
00226
00227       case DF_MATLAB_SPARSE:
00228         fprintf(file, "%% Size: %dx%d\n%% Nonzeros: %d\ntemp = zeros(%d, 3);\ntemp =[\n", this->size, this->size, Ap[this->size], Ap[this->size]);
00229         for (unsigned int j = 0; j < this->size; j++)
00230           for (unsigned int i = Ap[j]; i < Ap[j + 1]; i++)
00231           {
00232             fprintf(file, "%d %d ", Ai[i] + 1, j + 1);
00233             Hermes::Helpers::fprint_num(file, mumps_to_Scalar(Ax[i]));
00234             fprintf(file, "\n");
00235           }
00236           fprintf(file, "];\n%s = spconvert(temp);\n", var_name);
00237
00238           return true;
00239
00240       case DF_HERMES_BIN:
00241         {
00242           this->hermes_fwrite("HERMESX\001", 1, 8, file);
00243           int ssize = sizeof(Scalar);
00244           this->hermes_fwrite(&ssize, sizeof(int), 1, file);
00245           this->hermes_fwrite(&this->size, sizeof(int), 1, file);
00246           this->hermes_fwrite(&nnz, sizeof(int), 1, file);
00247           this->hermes_fwrite(Ap, sizeof(int), this->size + 1, file);
00248           this->hermes_fwrite(Ai, sizeof(int), nnz, file);
00249           this->hermes_fwrite(Ax, sizeof(Scalar), nnz, file);
00250           return true;
00251         }
00252
00253       default:
00254         return false;
00255       }
00256     }
00257
00258     template<typename Scalar>
00259     unsigned int MumpsMatrix<Scalar>::get_matrix_size() const
00260     {
00261       return this->size;
00262     }
00263
00264     template<typename Scalar>
00265     unsigned int MumpsMatrix<Scalar>::get_nnz() const
00266     {
00267       return nnz;
00268     }
00269
00270     template<typename Scalar>
00271     double MumpsMatrix<Scalar>::get_fill_in() const
00272     {
00273       return Ap[this->size] / (double) (this->size * this->size);
00274     }
00275
00276     template<typename Scalar>
00278     {
00280     };
00281
00282     template<typename Scalar>
00283     void MumpsMatrix<Scalar>::add_to_diagonal_blocks(int num_stages, MumpsMatrix<Scalar>* mat)
00284     {
00285       int ndof = mat->get_size();
00286       if(this->get_size() != (unsigned int) num_stages * ndof)
00287         throw Hermes::Exceptions::Exception("Incompatible matrix sizes in PetscMatrix<Scalar>::add_to_diagonal_blocks()");
00288
00289       for (int i = 0; i < num_stages; i++)
00290       {
00292       }
00293     }
00294
00295     template<typename Scalar>
00296     void MumpsMatrix<Scalar>::add_sparse_to_diagonal_blocks(int num_stages, SparseMatrix<Scalar>* mat)
00297     {
00299     }
00300
00301     inline ZMUMPS_COMPLEX& operator +=(ZMUMPS_COMPLEX &a, ZMUMPS_COMPLEX b)
00302     {
00303       a.r +=b.r;
00304       a.i +=b.i;
00305       return a;
00306     }
00307
00308     template<typename Scalar>
00309     void MumpsMatrix<Scalar>::add_as_block(unsigned int i, unsigned int j, MumpsMatrix<Scalar>* mat)
00310     {
00311       int idx;
00312       for (unsigned int col = 0;col<mat->get_size();col++)
00313       {
00314         for (unsigned int n = mat->Ap[col];n<mat->Ap[col + 1];n++)
00315         {
00316           idx = find_position(Ai + Ap[col + j], Ap[col + 1 + j] - Ap[col + j], mat->Ai[n] + i);
00317           if(idx<0)
00319           idx +=Ap[col + j];
00320           Ax[idx] +=mat->Ax[n];
00321         }
00322       }
00323     }
00324
00325     // Applies the matrix to vector_in and saves result to vector_out.
00326     template<typename Scalar>
00327     void MumpsMatrix<Scalar>::multiply_with_vector(Scalar* vector_in, Scalar* vector_out)
00328     {
00329       for(unsigned int i = 0;i<this->size;i++)
00330       {
00331         vector_out[i] = 0;
00332       }
00333       Scalar a;
00334       for (unsigned int i = 0;i<nnz;i++)
00335       {
00336         a = mumps_to_Scalar(Ax[i]);
00337         vector_out[jcn[i]-1] +=vector_in[irn[i]-1]*a;
00338       }
00339     }
00340     // Multiplies matrix with a Scalar.
00341     template<>
00342     void MumpsMatrix<double>::multiply_with_Scalar(double value)
00343     {
00344       int n = nnz;
00345       for(int i = 0;i<n;i++)
00346       {
00347         Ax[i] = Ax[i]*value;
00348       }
00349     }
00350
00351     template<>
00352     void MumpsMatrix<std::complex<double> >::multiply_with_Scalar(std::complex<double> value)
00353     {
00354       int n = nnz;
00355       std::complex<double> a;
00356       for(int i = 0;i<n;i++)
00357       {
00358         a = std::complex<double>(Ax[i].r, Ax[i].i);
00359         a = a*value;
00360         Ax[i].r = a.real();
00361         Ax[i].i = a.imag();
00362       }
00363     }
00364
00365     inline void mumps_assign_Scalar(ZMUMPS_COMPLEX & a, std::complex<double> b)
00366     {
00367       a.r = b.real();
00368       a.i = b.imag();
00369     }
00370
00371     inline void mumps_assign_Scalar(double & a, double b)
00372     {
00373       a = b;
00374     }
00375
00376     // Creates matrix using size, nnz, and the three arrays.
00377     template<typename Scalar>
00378     void MumpsMatrix<Scalar>::create(unsigned int size, unsigned int nnz, int* ap, int* ai, Scalar* ax)
00379     {
00380       this->nnz = nnz;
00381       this->size = size;
00382       this->Ap = new unsigned int[this->size + 1]; assert(this->Ap != NULL);
00383       this->Ai = new int[nnz];    assert(this->Ai != NULL);
00384       this->Ax = new typename mumps_type<Scalar>::mumps_Scalar[nnz]; assert(this->Ax != NULL);
00385       irn = new int[nnz];           assert(this->irn !=NULL);     // Row indices.
00386       jcn = new int[nnz];           assert(this->jcn !=NULL);     // Column indices.
00387
00388       for (unsigned int i = 0; i < this->size; i++)
00389       {
00390         this->Ap[i] = ap[i];
00391         for (int j = ap[i];j<ap[i + 1];j++) jcn[j] = i;
00392       }
00393       this->Ap[this->size] = ap[this->size];
00394       for (unsigned int i = 0; i < nnz; i++)
00395       {
00396         mumps_assign_Scalar(this->Ax[i], ax[i]);
00397         this->Ai[i] = ai[i];
00398         irn[i] = ai[i];
00399       }
00400     }
00401     // Duplicates a matrix (including allocation).
00402     template<typename Scalar>
00403     MumpsMatrix<Scalar>* MumpsMatrix<Scalar>::duplicate()
00404     {
00405       MumpsMatrix<Scalar> * nmat = new MumpsMatrix<Scalar>();
00406
00407       nmat->nnz = nnz;
00408       nmat->size = this->size;
00409       nmat->Ap = new unsigned int[this->size + 1]; assert(nmat->Ap != NULL);
00410       nmat->Ai = new int[nnz];    assert(nmat->Ai != NULL);
00411       nmat->Ax = new typename mumps_type<Scalar>::mumps_Scalar[nnz]; assert(nmat->Ax != NULL);
00412       nmat->irn = new int[nnz];           assert(nmat->irn !=NULL);     // Row indices.
00413       nmat->jcn = new int[nnz];           assert(nmat->jcn !=NULL);     // Column indices.
00414       for (unsigned int i = 0;i<nnz;i++)
00415       {
00416         nmat->Ai[i] = Ai[i];
00417         nmat->Ax[i] = Ax[i];
00418         nmat->irn[i] = irn[i];
00419         nmat->jcn[i] = jcn[i];
00420       }
00421       for (unsigned int i = 0;i<this->size + 1;i++)
00422       {
00423         nmat->Ap[i] = Ap[i];
00424       }
00425       return nmat;
00426     }
00427
00428     // MumpsVector<Scalar> /////////////////////////////////////////////////////////////////////////////////////
00429
00430     template<typename Scalar>
00431     MumpsVector<Scalar>::MumpsVector()
00432     {
00433       v = NULL;
00434       this->size = 0;
00435     }
00436
00437     template<typename Scalar>
00438     MumpsVector<Scalar>::~MumpsVector()
00439     {
00440       free();
00441     }
00442
00443     template<typename Scalar>
00444     void MumpsVector<Scalar>::alloc(unsigned int n)
00445     {
00446       free();
00447       this->size = n;
00448       v = new Scalar[n];
00449       zero();
00450     }
00451
00452     template<typename Scalar>
00453     void MumpsVector<Scalar>::change_sign()
00454     {
00455       for (unsigned int i = 0; i < this->size; i++) v[i] *= -1.;
00456     }
00457
00458     template<typename Scalar>
00459     void MumpsVector<Scalar>::zero()
00460     {
00461       memset(v, 0, this->size * sizeof(Scalar));
00462     }
00463
00464     template<typename Scalar>
00465     void MumpsVector<Scalar>::free()
00466     {
00467       delete [] v;
00468       v = NULL;
00469       this->size = 0;
00470     }
00471
00472     template<typename Scalar>
00473     void MumpsVector<Scalar>::set(unsigned int idx, Scalar y)
00474     {
00475       v[idx] = y;
00476     }
00477
00478     template<typename Scalar>
00479     void MumpsVector<Scalar>::add(unsigned int idx, Scalar y)
00480     {
00482       v[idx] += y;
00483     }
00484
00485     template<typename Scalar>
00486     void MumpsVector<Scalar>::add(unsigned int n, unsigned int *idx, Scalar *y)
00487     {
00488       for (unsigned int i = 0; i < n; i++)
00489       {
00490         v[idx[i]] += y[i];
00491       }
00492     }
00493
00494     template<typename Scalar>
00495     Scalar MumpsVector<Scalar>::get(unsigned int idx)
00496     {
00497       return v[idx];
00498     }
00499
00500     template<typename Scalar>
00501     void MumpsVector<Scalar>::extract(Scalar *v) const
00502     {
00503       memcpy(v, this->v, this->size * sizeof(Scalar));
00504     }
00505
00506     template<typename Scalar>
00508     {
00509       assert(this->length() == vec->length());
00510       for (unsigned int i = 0; i < this->length(); i++) this->add(i, vec->get(i));
00511     }
00512
00513     template<typename Scalar>
00515     {
00516       for (unsigned int i = 0; i < this->length(); i++) this->add(i, vec[i]);
00517     }
00518
00519     template<typename Scalar>
00520     bool MumpsVector<Scalar>::dump(FILE *file, const char *var_name, EMatrixDumpFormat fmt, char* number_format)
00521     {
00522       switch (fmt)
00523       {
00524       case DF_PLAIN_ASCII:
00525         for (unsigned int i = 0; i < this->size; i++)
00526         {
00527           Hermes::Helpers::fprint_num(file, v[i], number_format);
00528           fprintf(file, "\n");
00529         }
00530
00531         return true;
00532
00533       case DF_MATLAB_SPARSE:
00534         fprintf(file, "%% Size: %dx1\n%s =[\n", this->size, var_name);
00535         for (unsigned int i = 0; i < this->size; i++)
00536         {
00537           Hermes::Helpers::fprint_num(file, v[i], number_format);
00538           fprintf(file, "\n");
00539         }
00540         fprintf(file, " ];\n");
00541         return true;
00542
00543       case DF_HERMES_BIN:
00544         {
00545           this->hermes_fwrite("HERMESR\001", 1, 8, file);
00546           int ssize = sizeof(Scalar);
00547           this->hermes_fwrite(&ssize, sizeof(int), 1, file);
00548           this->hermes_fwrite(&this->size, sizeof(int), 1, file);
00549           this->hermes_fwrite(v, sizeof(Scalar), this->size, file);
00550           return true;
00551         }
00552
00553       default:
00554         return false;
00555       }
00556     }
00557
00558     template class HERMES_API MumpsMatrix<double>;
00559     template class HERMES_API MumpsMatrix<std::complex<double> >;
00560     template class HERMES_API MumpsVector<double>;
00561     template class HERMES_API MumpsVector<std::complex<double> >;
00562   }
00563   namespace Solvers
00564   {
00566 #define ICNTL(I)            icntl[(I)-1]
00567 #define MUMPS_INFO(param, I) (param).infog[(I)-1]
00568 #define INFOG(I)            infog[(I)-1]
00569
00571 #define JOB_INIT                    -1
00572 #define JOB_END                     -2
00573 #define JOB_ANALYZE_FACTORIZE_SOLVE  6
00574 #define JOB_FACTORIZE_SOLVE          5
00575 #define JOB_SOLVE                    3
00576
00577     template<>
00578     void MumpsSolver<double>::mumps_c(mumps_type<double>::mumps_struct * param)
00579     {
00580       dmumps_c(param);
00581     }
00582
00583     template<>
00584     void MumpsSolver<std::complex<double> >::mumps_c(mumps_type<std::complex<double> >::mumps_struct * param)
00585     {
00586       zmumps_c(param);
00587     }
00588
00589     template<typename Scalar>
00590     bool MumpsSolver<Scalar>::check_status()
00591     {
00592       switch (param.INFOG(1))
00593       {
00594       case 0: return true; // no error
00595       case -1: this->warn("Error occured on processor %d", MUMPS_INFO(param, 2)); break;
00597       default: this->warn("INFOG(1) = %d", param.INFOG(1)); break;
00598       }
00599       return false;
00600     }
00601
00602     template<typename Scalar>
00603     bool MumpsSolver<Scalar>::reinit()
00604     {
00605       if(inited)
00606       {
00607         // If there is already an instance of MUMPS running,
00608         // terminate it.
00609         param.job = JOB_END;
00610         mumps_c(&param);
00611       }
00612
00613       param.job = JOB_INIT;
00614       param.par = 1; // host also performs calculations
00615       param.sym = 0; // 0 = unsymmetric
00616       param.comm_fortran = USE_COMM_WORLD;
00617
00618       mumps_c(&param);
00619       inited = check_status();
00620
00621       if(inited)
00622       {
00623         // No printings.
00624         param.ICNTL(1) = -1;
00625         param.ICNTL(2) = -1;
00626         param.ICNTL(3) = -1;
00627         param.ICNTL(4) = 0;
00628
00629         param.ICNTL(20) = 0; // centralized dense RHS
00630         param.ICNTL(21) = 0; // centralized dense solution
00631
00632         // Specify the matrix.
00633         param.n = m->size;
00634         param.nz = m->nnz;
00635         param.irn = m->irn;
00636         param.jcn = m->jcn;
00637         param.a = m->Ax;
00638       }
00639
00640       return inited;
00641     }
00642
00643     template<typename Scalar>
00644     MumpsSolver<Scalar>::MumpsSolver(MumpsMatrix<Scalar> *m, MumpsVector<Scalar> *rhs) :
00645     DirectSolver<Scalar>(), m(m), rhs(rhs)
00646     {
00647       inited = false;
00648
00649       // Initial values for some fields of the MUMPS_STRUC structure that may be accessed
00650       // before MUMPS has been initialized.
00651       param.rhs = NULL;
00652       param.INFOG(33) = -999; // see the case HERMES_REUSE_MATRIX_REORDERING_AND_SCALING
00653       // in setup_factorization()
00654     }
00655
00656     template<typename Scalar>
00657     MumpsSolver<Scalar>::~MumpsSolver()
00658     {
00659       // Terminate the current instance of MUMPS.
00660       if(inited)
00661       {
00662         param.job = JOB_END;
00663         mumps_c(&param);
00664       }
00665
00666       if(param.rhs != NULL) delete [] param.rhs;
00667     }
00668
00669     template<typename Scalar>
00670     int MumpsSolver<Scalar>::get_matrix_size()
00671     {
00672       return m->size;
00673     }
00674
00675     template<typename Scalar>
00676     bool MumpsSolver<Scalar>::solve()
00677     {
00678       bool ret = false;
00679       assert(m != NULL);
00680       assert(rhs != NULL);
00681
00682       this->tick();
00683
00684       // Prepare the MUMPS data structure with input for the solver driver
00685       // (according to the chosen factorization reuse strategy), as well as
00686       // the system matrix.
00687       if( !setup_factorization() )
00688       {
00689         throw Hermes::Exceptions::LinearMatrixSolverException("LU factorization could not be completed.");
00690       }
00691
00692       // Specify the right-hand side (will be replaced by the solution).
00693       param.rhs = new typename mumps_type<Scalar>::mumps_Scalar[m->size];
00694       memcpy(param.rhs, rhs->v, m->size * sizeof(Scalar));
00695
00696       // Do the jobs specified in setup_factorization().
00697       mumps_c(&param);
00698
00699       ret = check_status();
00700
00701       if(ret)
00702       {
00703         delete [] this->sln;
00704         this->sln = new Scalar[m->size];
00705         for (unsigned int i = 0; i < rhs->size; i++)
00706           this->sln[i] = mumps_to_Scalar(param.rhs[i]);
00707       }
00708
00709       this->tick();
00710       this->time = this->accumulated();
00711
00712       delete [] param.rhs;
00713       param.rhs = NULL;
00714
00715       return ret;
00716     }
00717
00718     template<typename Scalar>
00719     bool MumpsSolver<Scalar>::setup_factorization()
00720     {
00721       // When called for the first time, all three phases (analysis, factorization,
00722       // solution) must be performed.
00723       int eff_fact_scheme = this->factorization_scheme;
00724       if(!inited)
00725         if( this->factorization_scheme == HERMES_REUSE_MATRIX_REORDERING ||
00726           this->factorization_scheme == HERMES_REUSE_FACTORIZATION_COMPLETELY )
00727           eff_fact_scheme = HERMES_FACTORIZE_FROM_SCRATCH;
00728
00729       switch (eff_fact_scheme)
00730       {
00731       case HERMES_FACTORIZE_FROM_SCRATCH:
00732         // (Re)initialize new instance.
00733         reinit();
00734
00735         // Let MUMPS decide when and how to compute matrix reordering and scaling.
00736         param.ICNTL(6) = 7;
00737         param.ICNTL(8) = 77;
00738         param.job = JOB_ANALYZE_FACTORIZE_SOLVE;
00739
00740         break;
00741       case HERMES_REUSE_MATRIX_REORDERING:
00742         // Let MUMPS reuse results of the symbolic analysis and perform
00743         // scaling during each factorization (values 1-8 may be set here,
00744         // corresponding to different scaling algorithms during factorization;
00745         // see the MUMPS documentation for details).
00746         param.ICNTL(8) = 7;
00747         param.job = JOB_FACTORIZE_SOLVE;
00748
00749         break;
00750       case HERMES_REUSE_MATRIX_REORDERING_AND_SCALING:
00751         // Perform scaling along with reordering during the symbolic analysis phase
00752         // and then reuse it during subsequent factorizations. New instance of MUMPS
00753         // has to be created before the analysis phase.
00754         if(param.INFOG(33) != -2)
00755         {
00756           reinit();
00757           param.ICNTL(6) = 5;
00758           param.job = JOB_ANALYZE_FACTORIZE_SOLVE;
00759           // After analysis is done, INFOG(33) will be set to -2 by MUMPS.
00760         }
00761         else
00762         {
00763           param.job = JOB_FACTORIZE_SOLVE;
00764         }
00765         break;
00766       case HERMES_REUSE_FACTORIZATION_COMPLETELY:
00767         param.job = JOB_SOLVE;
00768         break;
00769       }
00770
00771       return true;
00772     }
00773
00774     template class HERMES_API MumpsSolver<double>;
00775     template class HERMES_API MumpsSolver<std::complex<double> >;
00776   }
00777 }
00778 #endif
```