Amesos Package Browser (Single Doxygen Collection) Development
Amesos_MC64.cpp
Go to the documentation of this file.
00001 #include "Amesos_ConfigDefs.h"
00002 #include "Amesos_MC64.h"
00003 #include "Epetra_Comm.h"
00004 #include "Epetra_RowMatrix.h"
00005 #include <map>
00006 
00007 extern "C" void F77_FUNC(mc64id, MC64ID)(int*);
00008 extern "C" void F77_FUNC(mc64ad, MC64AD)(int*, int*, int*, int*, int*, 
00009                                          double*, int*, int*, int*, int*, 
00010                                          int*, double*, int*, int*);
00011 
00012 // ===========================================================================
00013 Amesos_MC64::Amesos_MC64(const Epetra_RowMatrix& A, int JOB,
00014                          const bool StoreTranspose, const bool analyze) :
00015   A_(A)
00016 {
00017   if (A_.Comm().NumProc() != 1)
00018   {
00019     std::cerr << "Class Amesos_MC64 can be used with one processor only!" << std::endl;
00020     exit(EXIT_FAILURE);
00021   }
00022   F77_FUNC(mc64id, MC64ID)(ICNTL_);
00023 
00024   Compute(JOB, StoreTranspose, analyze);
00025 }
00026 
00027 // ===========================================================================
00028 int Amesos_MC64::Compute(int JOB, const bool StoreTranspose, const bool analyze)
00029 {
00030   // convert A_ into column-based format
00031 
00032   int MaxNumEntries = A_.MaxNumEntries();
00033   int N = A_.NumMyRows();
00034   int NE = A_.NumMyNonzeros();
00035   std::vector<int> IP;
00036   std::vector<int> IRN;
00037   std::vector<double> VAL;
00038 
00039   std::vector<int> Indices(MaxNumEntries);
00040   std::vector<double> Values(MaxNumEntries);
00041 
00042   if (StoreTranspose)
00043   {
00044     // cheapest way, just store the transpose of A and not A. This is because
00045     // we have easy access to rows and not to columns.
00046     IP.resize(N + 1); IP[0] = 1;
00047     IRN.resize(NE);
00048     VAL.resize(NE);
00049     int count = 0;
00050     
00051     for (int i = 0 ; i < N ; ++i)
00052     {
00053       int NumEntries = 0;
00054 
00055       A_.ExtractMyRowCopy(i, MaxNumEntries, NumEntries, &Values[0], 
00056                           &Indices[0]);
00057 
00058       IP[i + 1] = IP[i] + NumEntries;
00059 
00060       for (int j = 0 ; j < NumEntries ; ++j)
00061       {
00062         IRN[count] = Indices[j] + 1;
00063         VAL[count] = Values[j];
00064         ++count;
00065       }
00066     }
00067     assert(count == NE);
00068   }
00069   else
00070   {
00071     // this stores the matrix and not its transpose, but it is more memory 
00072     // demading. The ifdef'd out part is simple and fast, but very memory
00073     // demanding.
00074 
00075 #if 0
00076     IRN.resize(N * MaxNumEntries);
00077     VAL.resize(N * MaxNumEntries);
00078 
00079     std::vector<int> count(N);
00080     for (int i = 0 ; i < N ; ++i) count[i] = 0;
00081 
00082     for (int i = 0 ; i < N ; ++i)
00083     {
00084       int NumEntries = 0;
00085 
00086       A_.ExtractMyRowCopy(i, MaxNumEntries, NumEntries, &Values[0], 
00087                           &Indices[0]);
00088 
00089       for (int j = 0 ; j < NumEntries ; ++j)
00090       {
00091         int col = Indices[j];
00092         IRN[col * MaxNumEntries + count[col]] = i + 1;
00093         VAL[col * MaxNumEntries + count[col]] = Values[j];
00094         ++count[col];
00095       }
00096     }
00097 
00098     // now compact storage
00099     int k = 0;
00100     for (int col = 0 ; col < N ; ++col)
00101     {
00102       for (int row = 0 ; row < count[col] ; ++row)
00103       {
00104         IRN[k] = IRN[col * MaxNumEntries + row];
00105         VAL[k] = VAL[col * MaxNumEntries + row];
00106         ++k;
00107       }
00108     }
00109     assert (k == NE);
00110 
00111     IRN.resize(k);
00112     VAL.resize(k);
00113 
00114     IP.resize(N + 1);
00115     IP[0] = 1;
00116 
00117     for (int col = 0 ; col < N ; ++col)
00118       IP[col + 1] = IP[col] + count[col];
00119 #else
00120     std::vector<std::vector<int> > cols(N);
00121     std::vector<std::vector<double> > vals(N);
00122 
00123     for (int i = 1 ; i <= N ; ++i)
00124     {
00125       int NumEntries = 0;
00126 
00127       A_.ExtractMyRowCopy(i - 1, MaxNumEntries, NumEntries, &Values[0], 
00128                           &Indices[0]);
00129 
00130       for (int j = 0 ; j < NumEntries ; ++j)
00131       {
00132         cols[Indices[j]].push_back(i);
00133         vals[Indices[j]].push_back(Values[j]);
00134       }
00135     }
00136 
00137     IP.resize(N + 1); IP[0] = 1;
00138     IRN.resize(NE);
00139     VAL.resize(NE);
00140     int count = 0;
00141 
00142     for (int i = 0 ; i < N ; ++i)
00143     {
00144       IP[i + 1] = IP[i] + cols[i].size();
00145 
00146       for (int j = 0 ; j < cols[i].size() ; ++j)
00147       {
00148         IRN[count] = cols[i][j];
00149         VAL[count] = vals[i][j];
00150         ++count;
00151       }
00152     }
00153 #endif
00154   }
00155 
00156   int NUM;
00157   CPERM_.resize(N);
00158   int LIW = 10 * N + NE;
00159   std::vector<int> IW(LIW);
00160   int LDW = 3 * N + NE;
00161   DW_.resize(LDW);
00162 
00163   JOB = 5;
00164   F77_FUNC(mc64ad, MC64aD)(&JOB, &N, &NE, &IP[0], &IRN[0], &VAL[0], 
00165                            &NUM, &CPERM_[0], &LIW, &IW[0], &LDW, 
00166                            &DW_[0], ICNTL_, INFO_);
00167 
00168   if (analyze)
00169   {
00170     std::map<double, int> table;
00171     for (int col = 0 ; col < N ; ++col)
00172     {
00173       for (int j = IP[col] ; j < IP[col + 1] ; ++j)
00174       {
00175         int row = IRN[j - 1] - 1;
00176         int new_col = CPERM_[col];
00177         double new_val = VAL[j - 1] * exp(DW_[row] + DW_[col + N]);
00178         if (new_val < 0.0) new_val = -new_val;
00179         if (new_val > 0.1) table[0.1]++;
00180         else if (new_val > 0.01) table[0.01]++;
00181         else if (new_val > 0.001) table[0.001]++;
00182         else if (new_val > 0.0001) table[0.0001]++;
00183         else if (new_val > 0.00001) table[0.00001]++;
00184         else if (new_val > 0.000001) table[0.000001]++;
00185         else table[0.0]++;
00186       }
00187     }
00188 
00189     std::cout << "# elements (total)    = " << NE << std::endl;
00190     std::cout << "# elements > 0.1      = " << table[0.1] << std::endl;
00191     std::cout << "# elements > 0.01     = " << table[0.01] << std::endl;
00192     std::cout << "# elements > 0.001    = " << table[0.001] << std::endl;
00193     std::cout << "# elements > 0.0001   = " << table[0.0001] << std::endl;
00194     std::cout << "# elements > 0.00001  = " << table[0.00001] << std::endl;
00195     std::cout << "# elements > 0.000001 = " << table[0.000001] << std::endl;
00196     std::cout << "# elements <=0.000001 = " << table[0.0] << std::endl;
00197   }
00198 
00199   AMESOS_RETURN(INFO_[0]);
00200 }
00201 
00202 // ===========================================================================
00203 double* Amesos_MC64::GetColScaling()
00204 {
00205   return((double*)&DW_[0 + A_.NumMyRows()]);
00206 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines