New matrix generators

matgen/generate_matrix_ge.cc

@@ -24,6 +24,7 @@
 #include <cstdio>
 #include <cstdlib>
 #include <utility>
+#include <numeric>      // std::gcd
 
 namespace slate {
 
@@ -203,7 +204,7 @@ void generate_matrix(
             const scalar_t inner_const = pi/scalar_t(max_mn+1);
 
             entry_type orthog_entry = [ outer_const, inner_const ]( int64_t i, int64_t j) {
-                scalar_t a = scalar_t(i) * scalar_t(j) * inner_const;
+                scalar_t a = scalar_t(i+1) * scalar_t(j+1) * inner_const;
                 return outer_const * sin(a);
             };
             set( orthog_entry, A, opts );
@@ -228,7 +229,8 @@ void generate_matrix(
         case TestMatrixType::ris: {
             const int64_t max_mn = std::max(n, m);
             entry_type ris_entry = [max_mn]( int64_t i, int64_t j ) {
-                return 0.5 / ( max_mn - j - i + 1.5 );
+                // n-(j + 1)-(i + 1)+1.5 = n-j-i-0.5
+                return 0.5 / ( max_mn - j - i - 0.5);
             };
             set( ris_entry, A, opts );
             break;
@@ -289,6 +291,177 @@ void generate_matrix(
             generate_geevx( params, dist, cond, sigma_max, A, Sigma, seed, opts );
             break;
         }
+
+        case TestMatrixType::minij: {
+            entry_type minij_entry = []( int64_t i, int64_t j ) {
+                return std::min(i + 1, j + 1);
+            };
+            set( minij_entry, A, opts );
+            break;
+        }
+
+        case TestMatrixType::hilb: {
+            entry_type hilb_entry = []( int64_t i, int64_t j ) {
+                return 1.0 / (i + j + 1);
+            };
+            set( hilb_entry, A, opts );
+            break;
+        }
+
+         case TestMatrixType::frank: {
+            const int64_t max_mn = std::max(n, m);
+            entry_type frank_entry = [max_mn]( int64_t i, int64_t j ) {
+                if ((i - j) > 1) {
+                    return int64_t( 0 );
+                }
+                else if ((i - j) == 1) {
+                    return max_mn - j - 1;
+                }
+                else {
+                    return max_mn - j;
+                }
+            };
+            set( frank_entry, A, opts );
+            break;
+        }
+
+        case TestMatrixType::lehmer: {
+            entry_type lehmer_entry = []( int64_t i, int64_t j ) {
+                return double (std::min(i, j) + 1) / (std::max(i, j) + 1);
+            };
+            set( lehmer_entry, A, opts );
+            break;
+        }
+
+        case TestMatrixType::lotkin: {
+            entry_type lotkin_entry = []( int64_t i, int64_t j ) {
+                return (i == 0 ? 1.0 : (1.0 / (i + j + 1)));
+            };
+            set( lotkin_entry, A, opts );
+            break;
+        }
+
+        case TestMatrixType::redheff: {
+            entry_type redheff_entry = []( int64_t i, int64_t j ) {
+                return ((j+1) % (i+1) == 0 || j == 0 ? 1 : 0);
+            };
+            set( redheff_entry, A, opts );
+            break;
+        }
+
+        case TestMatrixType::triw: {
+            entry_type triw_entry = []( int64_t i, int64_t j ) {
+                if (i == j) {
+                    return 1;
+                }
+                else if (i > j) {
+                    return 0;
+                }
+                else {
+                    return -1;
+                }
+            };
+            set( triw_entry, A, opts );
+            break;
+        }
+
+        case TestMatrixType::pei: {
+            entry_type pei_entry = []( int64_t i, int64_t j ) {
+                return (i == j ? 2 : 1);
+            };
+            set( pei_entry, A, opts);
+            break;
+        }
+
+        case TestMatrixType::tridiag: {
+            entry_type tridiag_entry = []( int64_t i, int64_t j ) {
+                if (i == j) {
+                    return 2;
+                }
+                else if (std::abs(i - j) == 1) {
+                    return -1;
+                }
+                else {
+                    return 0;
+                }
+            };
+            set( tridiag_entry, A, opts );
+            break;
+        }
+
+        case TestMatrixType::toeppen: {
+            entry_type toeppen_entry = []( int64_t i, int64_t j ) {
+                if (std::abs(j - i) == 1) {
+                    return int ((j - i) * 10);
+                }
+                else if (std::abs(i - j) == 2) {
+                    return 1;
+                }
+                else {
+                    return 0;
+                }
+            };
+            set( toeppen_entry, A, opts );
+            break;
+        }
+
+        case TestMatrixType::parter: {
+            entry_type parter_entry = []( int64_t i, int64_t j ) {
+                return 1 / (i - j + 0.5);
+            };
+            set( parter_entry, A, opts );
+            break;
+        }
+
+        case TestMatrixType::moler: {
+            entry_type moler_entry = []( int64_t i, int64_t j ) {
+                return (i == j ? i + 1 : std::min(i, j) - 1);
+            };
+            set( moler_entry, A, opts );
+            break;
+        }
+
+        case TestMatrixType::cauchy: {
+            entry_type cauchy_entry = []( int64_t i, int64_t j ) {
+                return 1.0 / ( i + j + 2);
+            };
+            set( cauchy_entry, A, opts );
+            break;
+        }
+
+        case TestMatrixType::chow: {
+            entry_type chow_entry = []( int64_t i, int64_t j ) {
+                return ((i - j) < -1 ? 0 : 1);
+            };
+            set( chow_entry, A, opts );
+            break;
+        }
+
+        case TestMatrixType::clement: {
+            const int64_t max_mn = std::max(n, m);
+            entry_type clement_entry = [max_mn]( int64_t i, int64_t j ) {
+                if ((i - j) == 1) {
+                    return max_mn - j - 1;
+                }
+                else if ((i - j) == -1) {
+                    return j;
+                }
+                else {
+                    return (int64_t)0;
+                }
+            };
+            set( clement_entry, A, opts );
+            break;
+        }
+
+        case TestMatrixType::gcdmat: {
+            entry_type gcdmat_entry = []( int64_t i, int64_t j ) {
+                return std::gcd(i + 1, j + 1);
+            };
+            set( gcdmat_entry, A, opts );
+            break;
+        }
+
     }
 
     // rand types have already been made diagonally dominant.

matgen/generate_matrix_utils.cc

@@ -84,6 +84,22 @@ void generate_matrix_usage()
     "riemann   |  matrix entry i,j equal to i+1 if j+2 divides i+2 else -1\n"
     "ris       |  matrix entry i,j equal to 0.5/(n-i-j+1.5)\n"
     "zielkeNS  |  nonsymmetric matrix of Zielke\n"
+    "minij     |  matrix entry i,j equal to min(i,j)\n"
+    "hilb      |  matrix entry i,j equal to 1/(i+j-1)\n"
+    "frank     |  square upper Hessenberg matrix with determinant 1\n"
+    "lehmer    |  matrix entry i,j equal to i/j if j>=i else j/i \n"
+    "lotkin    |  Hilbert matrix with 1's in its first row"
+    "redheff   |  matrix entry i,j equal to 1 if j==1 or i divides j, else 0\n"
+    "triw      |  upper triangular matrix with 1's on the diagonal and -1's above the diagonal\n"
+    "pei       |  2's on diagonal, rest 1's\n"
+    "tridiag   |  2's on diagonal, -1's on sub- and superdiagonal, rest zero\n"
+    "toeppen   |  1's on second subdiagonal and second superdiagonal, -10's on subdiagonal, 10's on superdiagonal, rest zero\n"
+    "parter    |  matrix entry i,j equal to 1/(i-j+0.5)\n"
+    "moler     |  matrix entry i,j equal to i on diagonal, rest min(i,j)\n"
+    "cauchy    |  matrix entry i,j equal to 1/(i + j)\n"
+    "chow      |  matrix with 1's on the diagonal, superdiagonal, and all elements below them, and 0's above the superdiagonal\n"
+    "clement   |  matrix with decreasing values [n-1 : 1] on subdiagonal, and increasing values [1 : n-1] on superdiagonal\n"
+    "gcdmat    |  matrix entry i,j equal to gcd(i,j)\n"
     "          |  \n"
     "rand@     |  matrix entries random uniform on (0, 1)\n"
     "rands@    |  matrix entries random uniform on (-1, 1)\n"
@@ -230,6 +246,22 @@ void decode_matrix(
              base == "syev"    ) { type = TestMatrixType::heev;     }
     else if (base == "geevx"   ) { type = TestMatrixType::geevx;    }
     else if (base == "geev"    ) { type = TestMatrixType::geev;     }
+    else if (base == "minij"   ) { type = TestMatrixType::minij;    }
+    else if (base == "hilb"    ) { type = TestMatrixType::hilb;     }
+    else if (base == "frank"   ) { type = TestMatrixType::frank;    }
+    else if (base == "lehmer"  ) { type = TestMatrixType::lehmer;   }
+    else if (base == "lotkin"  ) { type = TestMatrixType::lotkin;   }
+    else if (base == "redheff" ) { type = TestMatrixType::redheff;  }
+    else if (base == "triw"    ) { type = TestMatrixType::triw;     }
+    else if (base == "tridiag" ) { type = TestMatrixType::tridiag;  }
+    else if (base == "toeppen" ) { type = TestMatrixType::toeppen;  }
+    else if (base == "pei"     ) { type = TestMatrixType::pei;      }
+    else if (base == "parter"  ) { type = TestMatrixType::parter;   }
+    else if (base == "moler"   ) { type = TestMatrixType::moler;    }
+    else if (base == "cauchy"  ) { type = TestMatrixType::cauchy;   }
+    else if (base == "chow"    ) { type = TestMatrixType::chow;     }
+    else if (base == "clement" ) { type = TestMatrixType::clement;  }
+    else if (base == "gcdmat"  ) { type = TestMatrixType::gcdmat;   }
     else {
         snprintf( msg, sizeof( msg ), "in '%s': unknown matrix '%s'",
                   kind.c_str(), base.c_str() );

matgen/generate_matrix_utils.hh

@@ -53,6 +53,22 @@ enum class TestMatrixType {
     heev,
     geev,
     geevx,
+    minij,
+    hilb,
+    frank,
+    lehmer,
+    lotkin,
+    redheff,
+    triw,
+    tridiag,
+    toeppen,
+    pei,
+    parter,
+    moler,
+    cauchy,
+    chow,
+    clement,
+    gcdmat,
 };
 
 enum class TestMatrixDist {

test/run_tests.py

@@ -669,6 +669,8 @@ def filter_csv( values, csv ):
 
     [ 'scale_row_col', gen + dtype + mn + equed + nonuniform_nb + ge_matrix ],
 
+    [ 'matgen',    gen + dtype + mn + ab + nonuniform_nb + ge_matrix        ],
+
     [ 'set',    gen + dtype + mn + ab + nonuniform_nb + ge_matrix        ],
     [ 'tzset',  gen + dtype + mn + ab + nonuniform_nb + ge_matrix + uplo ],
     [ 'trset',  gen + dtype +  n + ab + nonuniform_nb + ge_matrix + uplo ],

test/test.cc

@@ -322,6 +322,8 @@ std::vector< testsweeper::routines_t > routines = {
     { "scale_row_col",      test_scale_row_col, Section::aux },
     { "",                   nullptr,           Section::newline },
 
+    { "matgen",             test_matgen,       Section::aux_gen },
+
     { "set",                test_set,          Section::aux },
     { "tzset",              test_set,          Section::aux },
     { "trset",              test_set,          Section::aux },

test/test.hh

@@ -448,6 +448,7 @@ void test_copy   (Params& params, bool run);
 void test_scale  (Params& params, bool run);
 void test_scale_row_col(Params& params, bool run);
 void test_set    (Params& params, bool run);
+void test_matgen (Params& params, bool run);
 
 //------------------------------------------------------------------------------
 inline double barrier_get_wtime(MPI_Comm comm)