MatSolve ordering support

Examples/MatSolve/main.cpp

@@ -14,6 +14,148 @@ using namespace INMOST;
 #define BARRIER
 #endif
 
+// developed by Igor Konshin
+#define NMAX 1024 //1024 // maximal pixel field dimension
+
+static int f[NMAX][NMAX]; // pixel field
+static int nm; // actual matrix dimention
+static int nf; // actual field dimention
+static double sc; // scaling indices from matrix to pixel field
+
+// Initialize print of matrix portrait into PS file
+void ps_init (int n)
+{
+    nm = n;
+    nf = (n <= NMAX) ? n : NMAX;
+    sc = 1;
+    if (nm > 1 && nm > nf) sc = (double)(nf-1)/(double)(nm-1);
+    std::cout<<"::: nm="<<nm<<" nf="<<nf<<" sc="<<sc<<std::endl;//db!
+    for (int i=0; i<nf; i++)
+        for (int j=0; j<nf; j++)
+            f[i][j] = 0;
+}
+
+// Save 1 nonzero into some pixel (j and j are started from 0 like in C)
+inline void ps_ij (int i, int j)
+{
+    if (i >= 0 && i < nm && j >= 0 && j < nm)
+        f[(int)(j*sc)][nf-1-(int)(i*sc)] = 1;
+}
+
+// Write the collected matrix portrait into PS file
+void ps_file (std::string file, int nbl = 1, int * ibl = NULL)
+{
+    double s1 = 50.0;
+    double s = 500.0 / nf;
+
+    std::fstream output(file.c_str(),std::ios::out);
+    output << "%%BoundingBox: 0 0 600 600" << std::endl;
+    output << "/m {moveto} def % x y" << std::endl;
+    output << "/l {lineto} def % x y" << std::endl;
+    output << "/s {stroke} def % x y" << std::endl;
+    output << "/n {newpath} def % x y" << std::endl;
+    output << "/c {closepath} def % x y" << std::endl;
+    output << 0.1 << " setlinewidth" << std::endl;
+    output << "n " << s1+s*.5 << " " << s1+s*.5 << " m "
+           << s1+s*(nf+.5) << " " << s1+s*.5 << " l "
+           << s1+s*(nf+.5) << " " << s1+s*(nf+.5) << " l "
+           << s1+s*.5 << " " << s1+s*(nf+.5) << " l "
+           << s1+s*.5 << " " << s1+s*.5 << " l s c" << std::endl;
+    if (nbl > 1) {
+        for (int i=1; i<nbl; i++) {
+            double x = sc*ibl[i] + .5;
+            double y = nf - (sc*ibl[i] + .5) + 1.0;
+            output << "n " << s1+s*.5 << " " << s1+s*y << " m "
+                   << s1+s*(nf+.5) << " " << s1+s*y << " l s c" << std::endl;
+            output << "n " << s1+s*x << " " << s1+s*.5 << " m "
+                   << s1+s*x << " " << s1+s*(nf+.5) << " l s c" << std::endl;
+        }
+    }
+    double r = (log10(1.0*nf) + 1.0) / 4.0;
+    r = s * ((r < 1.0) ? r : 1.0) / 2.0;
+    r = (r > 1.0) ? r : 1.0;
+    //r = .2;
+    output << 2*r << " setlinewidth" << std::endl;
+    output << "/d {moveto currentpoint" << r << " 0 360 arc fill} def % x y" << std::endl;
+    output << "n" << std::endl;
+
+    std::stringstream ps(std::ios::in | std::ios::out);
+    //ps << std::scientific;
+    //ps.precision(2);
+
+    for (int i=0; i<nf; i++)
+        for (int j=0; j<nf; j++)
+            if (f[i][j] != 0)
+                ps << s1+s*(i+1)-r << " " << s1+s*(j+1) << " m "
+                   << s1+s*(i+1)+r << " " << s1+s*(j+1) << " l" << std::endl;
+
+    output << ps.rdbuf();
+    output << "s c" << std::endl;
+    output << "showpage" << std::endl;
+}
+
+// Print of matrix portrait into PS file
+void ps_crs (int n, int nbl, int * ibl, int * ia, int *ja, const std::string file)
+{
+    int ia0 = ia[0];
+    int ja0 = n;
+    for (int k=0; k<ia[nbl]-ia0; k++) ja0 = (ja0 < ja[k]) ? ja0 : ja[k];
+    ps_init (n);
+    for (int i=0; i<n; i++)
+        for (int k=ia[i]-ia0; k<ia[i+1]-ia0; k++)
+            ps_ij (i, ja[k]-ja0);
+    ps_file (file, nbl, ibl);
+}
+
+#define USE_INMOST
+#if defined(USE_INMOST)
+// Print of INMOST matrix portrait into PS file
+void ps_inmost (Sparse::Matrix & A, const std::string file)
+{
+    int n = A.Size();
+    int row = A.GetFirstIndex();
+    ps_init (n);
+    for (Sparse::Matrix::iterator it = A.Begin(); it != A.End(); ++it) {
+        for (Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
+            ps_ij (row,jt->first);
+        row++;
+    }
+    ps_file (file, NULL, NULL);
+}
+#endif
+
+//#define PS_MATR_SELFTEST_CRS
+#if defined (PS_MATR_SELFTEST_CRS)
+int main ()
+{
+    //       [  3 -1 -1  0 ]
+    //   A = [ -1  3  0 -1 ]
+    //       [ -1  0  3 -1 ]
+    //       [  0 -1 -1  3 ]
+    static int n=4; // nza=12;
+    int   ia[5]  = { 0, 3, 6, 9, 12 }; //ia[n+1]
+    int   ja[12] = { 1,2,3, 1,2,4, 1,3,4, 2,3,4 }; //ja[nza]
+    int nbl=2, ibl[3]={0,2,4};
+    ps_crs (n,nbl,ibl,ia,ja,"z.ps");
+}
+#endif
+
+//#define PS_MATR_SELFTEST_0
+#if defined (PS_MATR_SELFTEST_0)
+int main ()
+{
+    int nbl=3, ibl[4]={0,30,60,100};
+    ps_init (100);
+    ps_ij (0,0);
+    ps_ij (1,1);
+    ps_ij (2,2);
+    ps_ij (12,17);
+    ps_ij (79,89);
+    ps_ij (99,99);
+    ps_file ("z.ps",nbl,ibl);
+}
+#endif
+
 int main(int argc, char **argv) {
     int processRank = 0, processorsCount = 1;
 
@@ -29,12 +171,16 @@ int main(int argc, char **argv) {
         std::string vectorXFileName = "";
         std::string parametersFileName = "";
         std::string solverName = "inner_ilu2";
+        std::string orderingFileName = "";
+
 
         bool matrixFound = false;
         bool vectorBFound = false;
         bool vectorXFound = false;
         bool parametersFound = false;
         bool typeFound = false;
+        bool saveVector = false;
+        bool orderingFound = false;
 
         //Parse argv parameters
         if (argc == 1) goto helpMessage;
@@ -51,8 +197,10 @@ int main(int argc, char **argv) {
                     std::cout << "Optional: " << std::endl;
                     std::cout << "-b, --bvector <RHS vector file name>" << std::endl;
                     std::cout << "-x, --xvector <X vector file name>" << std::endl;
-                    std::cout << "-p, --parameters <Solver parameters file name>" << std::endl;
+                    std::cout << "-d, --database <Solver parameters file name>" << std::endl;
                     std::cout << "-t, --type <Solver type name>" << std::endl;
+                    std::cout << "-ord, --ordering <Ordering file name>" << std::endl;
+                    std::cout << "-s, --save" << std::endl;
                     std::cout << "  Available solvers:" << std::endl;
                     Solver::Initialize(NULL, NULL, NULL);
                     std::vector<std::string> availableSolvers = Solver::getAvailableSolvers();
@@ -102,8 +250,8 @@ int main(int argc, char **argv) {
                 i++;
                 continue;
             }
-            //Parameters file name found with -p or --parameters options
-            if (strcmp(argv[i], "-p") == 0 || strcmp(argv[i], "--parameters") == 0) {
+            //Parameters file name found with -d or --database options
+            if (strcmp(argv[i], "-d") == 0 || strcmp(argv[i], "--database") == 0) {
                 if (processRank == 0) {
                     std::cout << "Solver parameters file found: " << argv[i + 1] << std::endl;
                 }
@@ -112,6 +260,20 @@ int main(int argc, char **argv) {
                 i++;
                 continue;
             }
+            //Ordering file name found with -ord or --parameters options
+            if (strcmp(argv[i], "-ord") == 0 || strcmp(argv[i], "--ordering") == 0) {
+                if (processRank == 0) {
+                    std::cout << "Ordering file found: " << argv[i + 1] << std::endl;
+                }
+                orderingFound = true;
+                orderingFileName = std::string(argv[i + 1]);
+                i++;
+                continue;
+            }
+            if (strcmp(argv[i], "-s") == 0 || strcmp(argv[i], "--save") == 0) {
+                saveVector = true;
+                continue;
+            }
             //Solver type found with -t ot --type options
             if (strcmp(argv[i], "-t") == 0 || strcmp(argv[i], "--type") == 0) {
                 if (processRank == 0) {
@@ -171,14 +333,28 @@ int main(int argc, char **argv) {
         Sparse::Vector b("rhs"); // Declare the right-hand side vector
         Sparse::Vector x("sol"); // Declare the solution vector
         double tempTimer = Timer(), solvingTimer;
+        if (orderingFound) {
+            if (processRank == 0) {
+                std::cout << "Using ordering file " << orderingFileName << std::endl;
+            }
+            mat.Load(matrixFileName, ENUMUNDEF, ENUMUNDEF, orderingFileName);
+        } else {
             mat.Load(matrixFileName); //if interval parameters not set, matrix will be divided automatically
+        }
+        if (processorsCount == 1) {
+            ps_inmost(mat, "a.ps"); //db !IK!
+        }
         BARRIER
         if (processRank == 0) {
             std::cout << "load matrix:    " << Timer() - tempTimer << std::endl;
         }
         tempTimer = Timer();
         if (vectorBFound) {
-            b.Load(vectorBFileName); // Load RHS vector
+            if (orderingFound) {
+                b.Load(vectorBFileName, ENUMUNDEF, ENUMUNDEF, orderingFileName); // Load RHS vector with ordering
+            } else {
+                b.Load(vectorBFileName); //if interval parameters not set, matrix will be divided automatically
+            }
             //b.Save(vectorBFileName + "_saved_test.rhs");
         } else { // Set local RHS to 1 if it was not specified
             INMOST_DATA_ENUM_TYPE mbeg, mend, k;
@@ -210,7 +386,10 @@ int main(int argc, char **argv) {
         iters = solver.Iterations(); // Get the number of iterations performed
         resid = solver.Residual();   // Get the final residual achieved
         reason = solver.ReturnReason(); // Get the convergence reason
-        //x.Save("output.sol");  // Save the solution if required
+
+        if (saveVector) {
+            x.Save("output.sol");  // Save the solution if required
+        }
 
         // Compute the true residual
         double aresid = 0, bresid = 0;

Examples/MatSolve/ps_matr.cpp

+// developed by Igor Konshin
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <cmath>
+
+#define NMAX 1024 // maximal pixel field dimention
+
+static int f[NMAX][NMAX]; // pixel field
+static int nm; // actual matrix dimention
+static int nf; // actual field dimention
+static double sc; // scaling indices from matrix to pixel field
+
+// Initialize print of matrix portrait into PS file
+void ps_init (int n)
+{
+    nm = n;
+    nf = (n <= NMAX) ? n : NMAX;
+    sc = 1;
+    if (nm > 1 && nm > nf) sc = (double)(nf-1)/(double)(nm-1);
+    std::cout<<"::: nm="<<nm<<" nf="<<nf<<" sc="<<sc<<std::endl;//db!
+    for (int i=0; i<nf; i++)
+        for (int j=0; j<nf; j++)
+            f[i][j] = 0;
+}
+
+// Save 1 nonzero into some pixel (j and j are started from 0 like in C)
+inline void ps_ij (int i, int j)
+{
+    if (i >= 0 && i < nm && j >= 0 && j < nm)
+        f[(int)(j*sc)][nf-1-(int)(i*sc)] = 1;
+}
+
+// Write the collected matrix portrait into PS file
+void ps_file (std::string file, int nbl = 1, int * ibl = NULL)
+{
+    double s1 = 50.0;
+    double s = 500.0 / nf;
+
+    std::fstream output(file.c_str(),std::ios::out);
+    output << "%%BoundingBox: 0 0 600 600" << std::endl;
+    output << "/m {moveto} def % x y" << std::endl;
+    output << "/l {lineto} def % x y" << std::endl;
+    output << "/s {stroke} def % x y" << std::endl;
+    output << "/n {newpath} def % x y" << std::endl;
+    output << "/c {closepath} def % x y" << std::endl;
+    output << 0.1 << " setlinewidth" << std::endl;
+    output << "n " << s1+s*.5 << " " << s1+s*.5 << " m "
+           << s1+s*(nf+.5) << " " << s1+s*.5 << " l "
+           << s1+s*(nf+.5) << " " << s1+s*(nf+.5) << " l "
+           << s1+s*.5 << " " << s1+s*(nf+.5) << " l "
+           << s1+s*.5 << " " << s1+s*.5 << " l s c" << std::endl;
+    if (nbl > 1) {
+        for (int i=1; i<nbl; i++) {
+            double x = sc*ibl[i] + .5;
+            double y = nf - (sc*ibl[i] + .5) + 1.0;
+            output << "n " << s1+s*.5 << " " << s1+s*y << " m "
+                   << s1+s*(nf+.5) << " " << s1+s*y << " l s c" << std::endl;
+            output << "n " << s1+s*x << " " << s1+s*.5 << " m "
+                   << s1+s*x << " " << s1+s*(nf+.5) << " l s c" << std::endl;
+        }
+    }
+    double r = (log10(1.0*nf) + 1.0) / 4.0;
+    r = s * ((r < 1.0) ? r : 1.0) / 2.0;
+    r = (r > 1.0) ? r : 1.0;
+    //r = .2;
+    output << 2*r << " setlinewidth" << std::endl;
+    output << "/d {moveto currentpoint" << r << " 0 360 arc fill} def % x y" << std::endl;
+    output << "n" << std::endl;
+
+    std::stringstream ps(std::ios::in | std::ios::out);
+    //ps << std::scientific;
+    //ps.precision(2);
+
+    for (int i=0; i<nf; i++)
+        for (int j=0; j<nf; j++)
+            if (f[i][j] != 0)
+                ps << s1+s*(i+1)-r << " " << s1+s*(j+1) << " m "
+                   << s1+s*(i+1)+r << " " << s1+s*(j+1) << " l" << std::endl;
+
+    output << ps.rdbuf();
+    output << "s c" << std::endl;
+    output << "showpage" << std::endl;
+}
+
+// Print of matrix portrait into PS file
+void ps_crs (int n, int nbl, int * ibl, int * ia, int *ja, const std::string file)
+{
+    int ia0 = ia[0];
+    int ja0 = n;
+    for (int k=0; k<ia[nbl]-ia0; k++) ja0 = (ja0 < ja[k]) ? ja0 : ja[k];
+    ps_init (n);
+    for (int i=0; i<n; i++)
+        for (int k=ia[i]-ia0; k<ia[i+1]-ia0; k++)
+            ps_ij (i, ja[k]-ja0);
+    ps_file (file, nbl, ibl);
+}
+
+#define USE_INMOST
+#if defined(USE_INMOST)
+// Print of INMOST matrix portrait into PS file
+void ps_inmost (Sparse::Matrix & A, const std::string file)
+{
+    int n = A.Size();
+    int row = A.GetFirstIndex();
+    ps_init (n);
+    for (Sparse::Matrix::iterator it = A.Begin(); it != A.End(); ++it) {
+        for (Solver::Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
+            ps_ij (row,jt->first);
+        row++;
+    }
+    ps_file (file, NULL, NULL);
+}
+#endif
+
+//#define PS_MATR_SELFTEST_CRS
+#if defined (PS_MATR_SELFTEST_CRS)
+int main ()
+{
+    //       [  3 -1 -1  0 ]
+    //   A = [ -1  3  0 -1 ]
+    //       [ -1  0  3 -1 ]
+    //       [  0 -1 -1  3 ]
+    static int n=4; // nza=12;
+    int   ia[5]  = { 0, 3, 6, 9, 12 }; //ia[n+1]
+    int   ja[12] = { 1,2,3, 1,2,4, 1,3,4, 2,3,4 }; //ja[nza]
+    int nbl=2, ibl[3]={0,2,4};
+    ps_crs (n,nbl,ibl,ia,ja,"z.ps");
+}
+#endif
+
+//#define PS_MATR_SELFTEST_0
+#if defined (PS_MATR_SELFTEST_0)
+int main ()
+{
+    int nbl=3, ibl[4]={0,30,60,100};
+    ps_init (100);
+    ps_ij (0,0);
+    ps_ij (1,1);
+    ps_ij (2,2);
+    ps_ij (12,17);
+    ps_ij (79,89);
+    ps_ij (99,99);
+    ps_file ("z.ps",nbl,ibl);
+}
+#endif

Source/Headers/inmost_sparse.h

@@ -100,7 +100,7 @@ namespace INMOST
             /// Load the vector from a single data file using the specified interval.
             /// If interval is not specified, then it will be automatically constructed,
             /// with the about equal block size (the last block may has larger dimension).
-            void                 Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg = ENUMUNDEF, INMOST_DATA_ENUM_TYPE mend = ENUMUNDEF);
+            void                 Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg = ENUMUNDEF, INMOST_DATA_ENUM_TYPE mend = ENUMUNDEF, std::string file_ord = "");
             
             bool                 & isParallel() {return is_parallel;}
             /// Get the vector name specified in the main constructor.
@@ -417,7 +417,7 @@ namespace INMOST
             /// Load the matrix from a single data file in MTX format using the specified interval.
             /// If interval is not specified, then it will be automatically constructed,
             /// with the about equal block size (the last block may has larger dimension).
-            void				 Load(std::string file, INMOST_DATA_ENUM_TYPE beg = ENUMUNDEF, INMOST_DATA_ENUM_TYPE end = ENUMUNDEF);
+            void				 Load(std::string file, INMOST_DATA_ENUM_TYPE beg = ENUMUNDEF, INMOST_DATA_ENUM_TYPE end = ENUMUNDEF, std::string file_ord = "");
             /// Save the distributed matrix to a single data file in MTX format using parallel MPI I/O.
             /// @see http://math.nist.gov/MatrixMarket/formats.html
             void                 Save(std::string file, const AnnotationService * annotation = NULL);

Source/Solver/sparse.cpp

@@ -662,7 +662,7 @@ namespace INMOST
         {
         }
 
-        void     Vector::Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg, INMOST_DATA_ENUM_TYPE mend)
+        void     Vector::Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg, INMOST_DATA_ENUM_TYPE mend, std::string file_ord)
         {
             char str[16384];
             std::ifstream input(file.c_str());
@@ -680,6 +680,30 @@ namespace INMOST
                 MPI_Comm_size(GetCommunicator(),&size);
             }
 #endif
+            int * ord = NULL;
+            if (file_ord != "")
+            {
+                std::ifstream input_ord;
+                input_ord.open(file_ord.c_str(), std::ifstream::in);
+                if( input_ord.fail() ) throw -2;
+                int n;
+                input_ord >> n;
+                ord = (int *) malloc(sizeof(int) * n);
+                for (int i=0; i<n; i++) input_ord >> ord[i];
+                int nbl;
+                input_ord >> nbl;
+                if( nbl != size ) throw -3;
+                int * ibl;
+                ibl = (int *) malloc(sizeof(int) * (nbl+1));
+                for (int i=0; i<nbl+1; i++) input_ord >> ibl[i];
+                if( mbeg == ENUMUNDEF ) mbeg = ibl[rank];
+                if( mend == ENUMUNDEF ) mend = ibl[rank+1];
+                for (int i=0; i<n; i++) ord[i] -= 1;
+                free(ibl);
+                input_ord.close();
+                //std::cout<<"Vector::Load(): n="<<n<<" nbl="<<nbl<<" np="<<size<<" id="<<rank<<" mbeg="<<mbeg<<" mend="<<mend<<std::endl;//db
+            }
+
             while( !input.getline(str,16384).eof() )
             {
                 k = 0; while( isspace(str[k]) ) k++;
@@ -700,12 +724,17 @@ namespace INMOST
                         break;
                     case 1:
                         istr >> val;
+                        if( ord ) {
+                            if( ord[ind] >= mbeg && ord[ind] < mend ) data[ord[ind]] = val;
+                        } else {
                             if( ind >= mbeg && ind < mend ) data[ind] = val;
+                        }
                         ind++;
                         break;
                 }
             }
             input.close();
+            if (file_ord != "") free(ord);
         }
 
 
@@ -958,6 +987,7 @@ namespace INMOST
             for(iterator it = Begin(); it != End(); ++it) nonzero += it->Size();
             //~ std::cout << rank << " total nonzero " << max_lines << " my nonzero " << nonzero << std::endl;
             input.close();
+            if (file_ord != "") free(ord);
         }
 
         void     Matrix::Save(std::string file, const AnnotationService * text)
@@ -1149,7 +1179,7 @@ namespace INMOST
 
 ////////class HessianMatrix
 
-        void     HessianMatrix::Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg, INMOST_DATA_ENUM_TYPE mend)
+        void     HessianMatrix::Load(std::string file, INMOST_DATA_ENUM_TYPE mbeg, INMOST_DATA_ENUM_TYPE mend, std::string file_ord)
         {
             char str[16384];
             std::ifstream input(file.c_str());

Tests/solver_test000/CMakeLists.txt

@@ -71,6 +71,10 @@ if(USE_SOLVER_ANI)
 add_test(NAME solver_test000_serial_ani COMMAND $<TARGET_FILE:solver_test000> 0 ani)
 endif()
 
+if(SOLVER_DEFINITIONS MATCHES "^.*HAVE_SOLVER_FCBIILU2.*$")
+  add_test(NAME solver_test000_serial_fcbiilu2 COMMAND $<TARGET_FILE:solver_test000> 0 fcbiilu2)
+endif()
+
 if( USE_MPI )
 if( EXISTS ${MPIEXEC} )
 add_test(NAME solver_test000_parallel_normal_inner_ilu2        COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 0 inner_ilu2)
@@ -114,5 +118,10 @@ add_test(NAME solver_test000_parallel_permute2_trilinos_ml     COMMAND ${MPIEXEC
 #add_test(NAME solver_test000_parallel_permute1_trilinos_belos  COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 1 6)
 #add_test(NAME solver_test000_parallel_permute2_trilinos_belos  COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 2 6)