htswanalysis/MACS/lib/gsl/gsl-1.11/linalg/test.c

   1 /* linalg/test.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2005, 2006, 2007 Gerard Jungman, Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman
  21  */
  22 #include <config.h>
  23 #include <stdlib.h>
  24 #include <gsl/gsl_test.h>
  25 #include <gsl/gsl_math.h>
  26 #include <gsl/gsl_ieee_utils.h>
  27 #include <gsl/gsl_permute_vector.h>
  28 #include <gsl/gsl_blas.h>
  29 #include <gsl/gsl_complex_math.h>
  30 #include <gsl/gsl_linalg.h>
  31
  32 #define TEST_SVD_4X4 1
  33
  34 int check (double x, double actual, double eps);
  35 gsl_matrix * create_hilbert_matrix(unsigned long size);
  36 gsl_matrix * create_general_matrix(unsigned long size1, unsigned long size2);
  37 gsl_matrix * create_vandermonde_matrix(unsigned long size);
  38 gsl_matrix * create_moler_matrix(unsigned long size);
  39 gsl_matrix * create_row_matrix(unsigned long size1, unsigned long size2);
  40 gsl_matrix * create_2x2_matrix(double a11, double a12, double a21, double a22);
  41 gsl_matrix * create_diagonal_matrix(double a[], unsigned long size);
  42
  43 int test_matmult(void);
  44 int test_matmult_mod(void);
  45 int test_LU_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  46 int test_LU_solve(void);
  47 int test_LUc_solve_dim(const gsl_matrix_complex * m, const double * actual, double eps);
  48 int test_LUc_solve(void);
  49 int test_QR_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  50 int test_QR_solve(void);
  51 int test_QR_QRsolve_dim(const gsl_matrix * m, const double * actual, double eps);
  52 int test_QR_QRsolve(void);
  53 int test_QR_lssolve_dim(const gsl_matrix * m, const double * actual, double eps);
  54 int test_QR_lssolve(void);
  55 int test_QR_decomp_dim(const gsl_matrix * m, double eps);
  56 int test_QR_decomp(void);
  57 int test_QRPT_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  58 int test_QRPT_solve(void);
  59 int test_QRPT_QRsolve_dim(const gsl_matrix * m, const double * actual, double eps);
  60 int test_QRPT_QRsolve(void);
  61 int test_QRPT_decomp_dim(const gsl_matrix * m, double eps);
  62 int test_QRPT_decomp(void);
  63 int test_QR_update_dim(const gsl_matrix * m, double eps);
  64 int test_QR_update(void);
  65
  66 int test_LQ_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  67 int test_LQ_solve(void);
  68 int test_LQ_LQsolve_dim(const gsl_matrix * m, const double * actual, double eps);
  69 int test_LQ_LQsolve(void);
  70 int test_LQ_lssolve_dim(const gsl_matrix * m, const double * actual, double eps);
  71 int test_LQ_lssolve(void);
  72 int test_LQ_decomp_dim(const gsl_matrix * m, double eps);
  73 int test_LQ_decomp(void);
  74 int test_PTLQ_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  75 int test_PTLQ_solve(void);
  76 int test_PTLQ_LQsolve_dim(const gsl_matrix * m, const double * actual, double eps);
  77 int test_PTLQ_LQsolve(void);
  78 int test_PTLQ_decomp_dim(const gsl_matrix * m, double eps);
  79 int test_PTLQ_decomp(void);
  80 int test_LQ_update_dim(const gsl_matrix * m, double eps);
  81 int test_LQ_update(void);
  82
  83 int test_SV_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  84 int test_SV_solve(void);
  85 int test_SV_decomp_dim(const gsl_matrix * m, double eps);
  86 int test_SV_decomp(void);
  87 int test_SV_decomp_mod_dim(const gsl_matrix * m, double eps);
  88 int test_SV_decomp_mod(void);
  89 int test_SV_decomp_jacobi_dim(const gsl_matrix * m, double eps);
  90 int test_SV_decomp_jacobi(void);
  91 int test_cholesky_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  92 int test_cholesky_solve(void);
  93 int test_cholesky_decomp_dim(const gsl_matrix * m, double eps);
  94 int test_cholesky_decomp(void);
  95 int test_HH_solve_dim(const gsl_matrix * m, const double * actual, double eps);
  96 int test_HH_solve(void);
  97 int test_TDS_solve_dim(unsigned long dim, double d, double od, const double * actual, double eps);
  98 int test_TDS_solve(void);
  99 int test_TDN_solve_dim(unsigned long dim, double d, double a, double b, const double * actual, double eps);
 100 int test_TDN_solve(void);
 101 int test_TDS_cyc_solve_one(const unsigned long dim, const double * d, const double * od, const double * r,
 102                           const double * actual, double eps);
 103 int test_TDS_cyc_solve(void);
 104 int test_TDN_cyc_solve_dim(unsigned long dim, double d, double a, double b, const double * actual, double eps);
 105 int test_TDN_cyc_solve(void);
 106 int test_bidiag_decomp_dim(const gsl_matrix * m, double eps);
 107 int test_bidiag_decomp(void);
 108
 109 int
 110 check (double x, double actual, double eps)
 111 {
 112   if (x == actual)
 113     {
 114       return 0;
 115     }
 116   else if (actual == 0)
 117     {
 118       return fabs(x) > eps;
 119     }
 120   else
 121     {
 122       return (fabs(x - actual)/fabs(actual)) > eps;
 123     }
 124 }
 125
 126 gsl_matrix *
 127 create_hilbert_matrix(unsigned long size)
 128 {
 129   unsigned long i, j;
 130   gsl_matrix * m = gsl_matrix_alloc(size, size);
 131   for(i=0; i<size; i++) {
 132     for(j=0; j<size; j++) {
 133       gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
 134     }
 135   }
 136   return m;
 137 }
 138
 139 gsl_matrix *
 140 create_general_matrix(unsigned long size1, unsigned long size2)
 141 {
 142   unsigned long i, j;
 143   gsl_matrix * m = gsl_matrix_alloc(size1, size2);
 144   for(i=0; i<size1; i++) {
 145     for(j=0; j<size2; j++) {
 146       gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
 147     }
 148   }
 149   return m;
 150 }
 151
 152 gsl_matrix *
 153 create_singular_matrix(unsigned long size1, unsigned long size2)
 154 {
 155   unsigned long i, j;
 156   gsl_matrix * m = gsl_matrix_alloc(size1, size2);
 157   for(i=0; i<size1; i++) {
 158     for(j=0; j<size2; j++) {
 159       gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
 160     }
 161   }
 162
 163   /* zero the first column */
 164   for(j = 0; j <size2; j++)
 165     gsl_matrix_set(m,0,j,0.0);
 166
 167   return m;
 168 }
 169
 170
 171 gsl_matrix *
 172 create_vandermonde_matrix(unsigned long size)
 173 {
 174   unsigned long i, j;
 175   gsl_matrix * m = gsl_matrix_alloc(size, size);
 176   for(i=0; i<size; i++) {
 177     for(j=0; j<size; j++) {
 178       gsl_matrix_set(m, i, j, pow(i + 1.0, size - j - 1.0));
 179     }
 180   }
 181   return m;
 182 }
 183
 184 gsl_matrix *
 185 create_moler_matrix(unsigned long size)
 186 {
 187   unsigned long i, j;
 188   gsl_matrix * m = gsl_matrix_alloc(size, size);
 189   for(i=0; i<size; i++) {
 190     for(j=0; j<size; j++) {
 191       gsl_matrix_set(m, i, j, GSL_MIN(i+1,j+1)-2.0);
 192     }
 193   }
 194   return m;
 195 }
 196
 197 gsl_matrix_complex *
 198 create_complex_matrix(unsigned long size)
 199 {
 200   unsigned long i, j;
 201   gsl_matrix_complex * m = gsl_matrix_complex_alloc(size, size);
 202   for(i=0; i<size; i++) {
 203     for(j=0; j<size; j++) {
 204       gsl_complex z = gsl_complex_rect(1.0/(i+j+1.0), 1/(i*i+j*j+0.5));
 205       gsl_matrix_complex_set(m, i, j, z);
 206     }
 207   }
 208   return m;
 209 }
 210
 211 gsl_matrix *
 212 create_row_matrix(unsigned long size1, unsigned long size2)
 213 {
 214   unsigned long i;
 215   gsl_matrix * m = gsl_matrix_calloc(size1, size2);
 216   for(i=0; i<size1; i++) {
 217       gsl_matrix_set(m, i, 0, 1.0/(i+1.0));
 218   }
 219
 220   return m;
 221 }
 222
 223 gsl_matrix *
 224 create_2x2_matrix(double a11, double a12, double a21, double a22)
 225 {
 226   gsl_matrix * m = gsl_matrix_alloc(2, 2);
 227   gsl_matrix_set(m, 0, 0, a11);
 228   gsl_matrix_set(m, 0, 1, a12);
 229   gsl_matrix_set(m, 1, 0, a21);
 230   gsl_matrix_set(m, 1, 1, a22);
 231   return m;
 232 }
 233
 234 gsl_matrix *
 235 create_diagonal_matrix(double a[], unsigned long size)
 236 {
 237   unsigned long i;
 238   gsl_matrix * m = gsl_matrix_calloc(size, size);
 239   for(i=0; i<size; i++) {
 240       gsl_matrix_set(m, i, i, a[i]);
 241   }
 242
 243   return m;
 244 }
 245
 246 gsl_matrix * m11;
 247 gsl_matrix * m51;
 248
 249 gsl_matrix * m35;
 250 gsl_matrix * m53;
 251 gsl_matrix * m97;
 252
 253 gsl_matrix * s35;
 254 gsl_matrix * s53;
 255
 256 gsl_matrix * hilb2;
 257 gsl_matrix * hilb3;
 258 gsl_matrix * hilb4;
 259 gsl_matrix * hilb12;
 260
 261 gsl_matrix * row3;
 262 gsl_matrix * row5;
 263 gsl_matrix * row12;
 264
 265 gsl_matrix * A22;
 266 gsl_matrix * A33;
 267 gsl_matrix * A44;
 268 gsl_matrix * A55;
 269
 270 gsl_matrix_complex * c7;
 271
 272 gsl_matrix * inf5; double inf5_data[] = {1.0, 0.0, -3.0, 0.0, -5.0};
 273 gsl_matrix * nan5;
 274 gsl_matrix * dblmin3, * dblmin5;
 275
 276 double m53_lssolution[] = {52.5992295702070, -337.7263113752073,
 277                            351.8823436427604};
 278 double hilb2_solution[] = {-8.0, 18.0} ;
 279 double hilb3_solution[] = {27.0, -192.0, 210.0};
 280 double hilb4_solution[] = {-64.0, 900.0, -2520.0, 1820.0};
 281 double hilb12_solution[] = {-1728.0, 245388.0, -8528520.0,
 282                             127026900.0, -1009008000.0, 4768571808.0,
 283                             -14202796608.0, 27336497760.0, -33921201600.0,
 284                             26189163000.0, -11437874448.0, 2157916488.0 };
 285
 286 double c7_solution[] = { 2.40717272023734e+01, -9.84612797621247e+00,
 287                          -2.69338853034031e+02, 8.75455232472528e+01,
 288                          2.96661356736296e+03, -1.02624473923993e+03,
 289                          -1.82073812124749e+04, 5.67384473042410e+03,
 290                          5.57693879019068e+04, -1.61540963210502e+04,
 291                          -7.88941207561151e+04, 1.95053812987858e+04,
 292                          3.95548551241728e+04, -7.76593696255317e+03 };
 293
 294 gsl_matrix * vander2;
 295 gsl_matrix * vander3;
 296 gsl_matrix * vander4;
 297 gsl_matrix * vander12;
 298
 299 double vander2_solution[] = {1.0, 0.0};
 300 double vander3_solution[] = {0.0, 1.0, 0.0};
 301 double vander4_solution[] = {0.0, 0.0, 1.0, 0.0};
 302 double vander12_solution[] = {0.0, 0.0, 0.0, 0.0,
 303                             0.0, 0.0, 0.0, 0.0,
 304                             0.0, 0.0, 1.0, 0.0};
 305
 306 gsl_matrix * moler10;
 307
 308 /* matmult now obsolete */
 309 #ifdef MATMULT
 310 int
 311 test_matmult(void)
 312 {
 313   int s = 0;
 314
 315   gsl_matrix * A = gsl_matrix_calloc(2, 2);
 316   gsl_matrix * B = gsl_matrix_calloc(2, 3);
 317   gsl_matrix * C = gsl_matrix_calloc(2, 3);
 318
 319   gsl_matrix_set(A, 0, 0, 10.0);
 320   gsl_matrix_set(A, 0, 1,  5.0);
 321   gsl_matrix_set(A, 1, 0,  1.0);
 322   gsl_matrix_set(A, 1, 1, 20.0);
 323
 324   gsl_matrix_set(B, 0, 0, 10.0);
 325   gsl_matrix_set(B, 0, 1,  5.0);
 326   gsl_matrix_set(B, 0, 2,  2.0);
 327   gsl_matrix_set(B, 1, 0,  1.0);
 328   gsl_matrix_set(B, 1, 1,  3.0);
 329   gsl_matrix_set(B, 1, 2,  2.0);
 330
 331   gsl_linalg_matmult(A, B, C);
 332
 333   s += ( fabs(gsl_matrix_get(C, 0, 0) - 105.0) > GSL_DBL_EPSILON );
 334   s += ( fabs(gsl_matrix_get(C, 0, 1) -  65.0) > GSL_DBL_EPSILON );
 335   s += ( fabs(gsl_matrix_get(C, 0, 2) -  30.0) > GSL_DBL_EPSILON );
 336   s += ( fabs(gsl_matrix_get(C, 1, 0) -  30.0) > GSL_DBL_EPSILON );
 337   s += ( fabs(gsl_matrix_get(C, 1, 1) -  65.0) > GSL_DBL_EPSILON );
 338   s += ( fabs(gsl_matrix_get(C, 1, 2) -  42.0) > GSL_DBL_EPSILON );
 339
 340   gsl_matrix_free(A);
 341   gsl_matrix_free(B);
 342   gsl_matrix_free(C);
 343
 344   return s;
 345 }
 346
 347
 348 int
 349 test_matmult_mod(void)
 350 {
 351   int s = 0;
 352
 353   gsl_matrix * A = gsl_matrix_calloc(3, 3);
 354   gsl_matrix * B = gsl_matrix_calloc(3, 3);
 355   gsl_matrix * C = gsl_matrix_calloc(3, 3);
 356   gsl_matrix * D = gsl_matrix_calloc(2, 3);
 357   gsl_matrix * E = gsl_matrix_calloc(2, 3);
 358   gsl_matrix * F = gsl_matrix_calloc(2, 2);
 359
 360   gsl_matrix_set(A, 0, 0, 10.0);
 361   gsl_matrix_set(A, 0, 1,  5.0);
 362   gsl_matrix_set(A, 0, 2,  1.0);
 363   gsl_matrix_set(A, 1, 0,  1.0);
 364   gsl_matrix_set(A, 1, 1, 20.0);
 365   gsl_matrix_set(A, 1, 2,  5.0);
 366   gsl_matrix_set(A, 2, 0,  1.0);
 367   gsl_matrix_set(A, 2, 1,  3.0);
 368   gsl_matrix_set(A, 2, 2,  7.0);
 369
 370   gsl_matrix_set(B, 0, 0, 10.0);
 371   gsl_matrix_set(B, 0, 1,  5.0);
 372   gsl_matrix_set(B, 0, 2,  2.0);
 373   gsl_matrix_set(B, 1, 0,  1.0);
 374   gsl_matrix_set(B, 1, 1,  3.0);
 375   gsl_matrix_set(B, 1, 2,  2.0);
 376   gsl_matrix_set(B, 2, 0,  1.0);
 377   gsl_matrix_set(B, 2, 1,  3.0);
 378   gsl_matrix_set(B, 2, 2,  2.0);
 379
 380   gsl_matrix_set(D, 0, 0, 10.0);
 381   gsl_matrix_set(D, 0, 1,  5.0);
 382   gsl_matrix_set(D, 0, 2,  1.0);
 383   gsl_matrix_set(D, 1, 0,  1.0);
 384   gsl_matrix_set(D, 1, 1, 20.0);
 385   gsl_matrix_set(D, 1, 2,  5.0);
 386
 387   gsl_matrix_set(E, 0, 0, 10.0);
 388   gsl_matrix_set(E, 0, 1,  5.0);
 389   gsl_matrix_set(E, 0, 2,  2.0);
 390   gsl_matrix_set(E, 1, 0,  1.0);
 391   gsl_matrix_set(E, 1, 1,  3.0);
 392   gsl_matrix_set(E, 1, 2,  2.0);
 393
 394   gsl_linalg_matmult_mod(A, GSL_LINALG_MOD_NONE, B, GSL_LINALG_MOD_NONE, C);
 395   s += ( fabs(gsl_matrix_get(C, 0, 0) - 106.0) > GSL_DBL_EPSILON );
 396   s += ( fabs(gsl_matrix_get(C, 0, 1) -  68.0) > GSL_DBL_EPSILON );
 397   s += ( fabs(gsl_matrix_get(C, 0, 2) -  32.0) > GSL_DBL_EPSILON );
 398   s += ( fabs(gsl_matrix_get(C, 1, 0) -  35.0) > GSL_DBL_EPSILON );
 399   s += ( fabs(gsl_matrix_get(C, 1, 1) -  80.0) > GSL_DBL_EPSILON );
 400   s += ( fabs(gsl_matrix_get(C, 1, 2) -  52.0) > GSL_DBL_EPSILON );
 401   s += ( fabs(gsl_matrix_get(C, 2, 0) -  20.0) > GSL_DBL_EPSILON );
 402   s += ( fabs(gsl_matrix_get(C, 2, 1) -  35.0) > GSL_DBL_EPSILON );
 403   s += ( fabs(gsl_matrix_get(C, 2, 2) -  22.0) > GSL_DBL_EPSILON );
 404
 405   gsl_linalg_matmult_mod(A, GSL_LINALG_MOD_TRANSPOSE, B, GSL_LINALG_MOD_NONE, C);
 406   s += ( fabs(gsl_matrix_get(C, 0, 0) - 102.0) > GSL_DBL_EPSILON );
 407   s += ( fabs(gsl_matrix_get(C, 0, 1) -  56.0) > GSL_DBL_EPSILON );
 408   s += ( fabs(gsl_matrix_get(C, 0, 2) -  24.0) > GSL_DBL_EPSILON );
 409   s += ( fabs(gsl_matrix_get(C, 1, 0) -  73.0) > GSL_DBL_EPSILON );
 410   s += ( fabs(gsl_matrix_get(C, 1, 1) -  94.0) > GSL_DBL_EPSILON );
 411   s += ( fabs(gsl_matrix_get(C, 1, 2) -  56.0) > GSL_DBL_EPSILON );
 412   s += ( fabs(gsl_matrix_get(C, 2, 0) -  22.0) > GSL_DBL_EPSILON );
 413   s += ( fabs(gsl_matrix_get(C, 2, 1) -  41.0) > GSL_DBL_EPSILON );
 414   s += ( fabs(gsl_matrix_get(C, 2, 2) -  26.0) > GSL_DBL_EPSILON );
 415
 416   gsl_linalg_matmult_mod(A, GSL_LINALG_MOD_NONE, B, GSL_LINALG_MOD_TRANSPOSE, C);
 417   s += ( fabs(gsl_matrix_get(C, 0, 0) - 127.0) > GSL_DBL_EPSILON );
 418   s += ( fabs(gsl_matrix_get(C, 0, 1) -  27.0) > GSL_DBL_EPSILON );
 419   s += ( fabs(gsl_matrix_get(C, 0, 2) -  27.0) > GSL_DBL_EPSILON );
 420   s += ( fabs(gsl_matrix_get(C, 1, 0) - 120.0) > GSL_DBL_EPSILON );
 421   s += ( fabs(gsl_matrix_get(C, 1, 1) -  71.0) > GSL_DBL_EPSILON );
 422   s += ( fabs(gsl_matrix_get(C, 1, 2) -  71.0) > GSL_DBL_EPSILON );
 423   s += ( fabs(gsl_matrix_get(C, 2, 0) -  39.0) > GSL_DBL_EPSILON );
 424   s += ( fabs(gsl_matrix_get(C, 2, 1) -  24.0) > GSL_DBL_EPSILON );
 425   s += ( fabs(gsl_matrix_get(C, 2, 2) -  24.0) > GSL_DBL_EPSILON );
 426
 427   gsl_linalg_matmult_mod(A, GSL_LINALG_MOD_TRANSPOSE, B, GSL_LINALG_MOD_TRANSPOSE, C);
 428   s += ( fabs(gsl_matrix_get(C, 0, 0) - 107.0) > GSL_DBL_EPSILON );
 429   s += ( fabs(gsl_matrix_get(C, 0, 1) -  15.0) > GSL_DBL_EPSILON );
 430   s += ( fabs(gsl_matrix_get(C, 0, 2) -  15.0) > GSL_DBL_EPSILON );
 431   s += ( fabs(gsl_matrix_get(C, 1, 0) - 156.0) > GSL_DBL_EPSILON );
 432   s += ( fabs(gsl_matrix_get(C, 1, 1) -  71.0) > GSL_DBL_EPSILON );
 433   s += ( fabs(gsl_matrix_get(C, 1, 2) -  71.0) > GSL_DBL_EPSILON );
 434   s += ( fabs(gsl_matrix_get(C, 2, 0) -  49.0) > GSL_DBL_EPSILON );
 435   s += ( fabs(gsl_matrix_get(C, 2, 1) -  30.0) > GSL_DBL_EPSILON );
 436   s += ( fabs(gsl_matrix_get(C, 2, 2) -  30.0) > GSL_DBL_EPSILON );
 437
 438   /* now try for non-symmetric matrices */
 439   gsl_linalg_matmult_mod(D, GSL_LINALG_MOD_TRANSPOSE, E, GSL_LINALG_MOD_NONE, C);
 440   s += ( fabs(gsl_matrix_get(C, 0, 0) - 101.0) > GSL_DBL_EPSILON );
 441   s += ( fabs(gsl_matrix_get(C, 0, 1) -  53.0) > GSL_DBL_EPSILON );
 442   s += ( fabs(gsl_matrix_get(C, 0, 2) -  22.0) > GSL_DBL_EPSILON );
 443   s += ( fabs(gsl_matrix_get(C, 1, 0) -  70.0) > GSL_DBL_EPSILON );
 444   s += ( fabs(gsl_matrix_get(C, 1, 1) -  85.0) > GSL_DBL_EPSILON );
 445   s += ( fabs(gsl_matrix_get(C, 1, 2) -  50.0) > GSL_DBL_EPSILON );
 446   s += ( fabs(gsl_matrix_get(C, 2, 0) -  15.0) > GSL_DBL_EPSILON );
 447   s += ( fabs(gsl_matrix_get(C, 2, 1) -  20.0) > GSL_DBL_EPSILON );
 448   s += ( fabs(gsl_matrix_get(C, 2, 2) -  12.0) > GSL_DBL_EPSILON );
 449
 450
 451   gsl_linalg_matmult_mod(D, GSL_LINALG_MOD_NONE, E, GSL_LINALG_MOD_TRANSPOSE, F);
 452   s += ( fabs(gsl_matrix_get(F, 0, 0) - 127.0) > GSL_DBL_EPSILON );
 453   s += ( fabs(gsl_matrix_get(F, 0, 1) -  27.0) > GSL_DBL_EPSILON );
 454   s += ( fabs(gsl_matrix_get(F, 1, 0) - 120.0) > GSL_DBL_EPSILON );
 455   s += ( fabs(gsl_matrix_get(F, 1, 1) -  71.0) > GSL_DBL_EPSILON );
 456
 457
 458   gsl_matrix_free(A);
 459   gsl_matrix_free(B);
 460   gsl_matrix_free(C);
 461   gsl_matrix_free(D);
 462   gsl_matrix_free(E);
 463   gsl_matrix_free(F);
 464
 465   return s;
 466 }
 467 #endif
 468
 469 int
 470 test_LU_solve_dim(const gsl_matrix * m, const double * actual, double eps)
 471 {
 472   int s = 0;
 473   int signum;
 474   unsigned long i, dim = m->size1;
 475
 476   gsl_permutation * perm = gsl_permutation_alloc(dim);
 477   gsl_vector * rhs = gsl_vector_alloc(dim);
 478   gsl_matrix * lu  = gsl_matrix_alloc(dim,dim);
 479   gsl_vector * x = gsl_vector_alloc(dim);
 480   gsl_vector * residual = gsl_vector_alloc(dim);
 481   gsl_matrix_memcpy(lu,m);
 482   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
 483   s += gsl_linalg_LU_decomp(lu, perm, &signum);
 484   s += gsl_linalg_LU_solve(lu, perm, rhs, x);
 485
 486   for(i=0; i<dim; i++) {
 487     int foo = check(gsl_vector_get(x, i),actual[i],eps);
 488     if(foo) {
 489       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
 490     }
 491     s += foo;
 492   }
 493
 494   s += gsl_linalg_LU_refine(m, lu, perm, rhs, x, residual);
 495
 496   for(i=0; i<dim; i++) {
 497     int foo = check(gsl_vector_get(x, i),actual[i],eps);
 498     if(foo) {
 499       printf("%3lu[%lu]: %22.18g   %22.18g (improved)\n", dim, i, gsl_vector_get(x, i), actual[i]);
 500     }
 501     s += foo;
 502   }
 503
 504   gsl_vector_free(residual);
 505   gsl_vector_free(x);
 506   gsl_matrix_free(lu);
 507   gsl_vector_free(rhs);
 508   gsl_permutation_free(perm);
 509
 510   return s;
 511 }
 512
 513
 514 int test_LU_solve(void)
 515 {
 516   int f;
 517   int s = 0;
 518
 519   f = test_LU_solve_dim(hilb2, hilb2_solution, 8.0 * GSL_DBL_EPSILON);
 520   gsl_test(f, "  LU_solve hilbert(2)");
 521   s += f;
 522
 523   f = test_LU_solve_dim(hilb3, hilb3_solution, 64.0 * GSL_DBL_EPSILON);
 524   gsl_test(f, "  LU_solve hilbert(3)");
 525   s += f;
 526
 527   f = test_LU_solve_dim(hilb4, hilb4_solution, 2048.0 * GSL_DBL_EPSILON);
 528   gsl_test(f, "  LU_solve hilbert(4)");
 529   s += f;
 530
 531   f = test_LU_solve_dim(hilb12, hilb12_solution, 0.5);
 532   gsl_test(f, "  LU_solve hilbert(12)");
 533   s += f;
 534
 535   f = test_LU_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
 536   gsl_test(f, "  LU_solve vander(2)");
 537   s += f;
 538
 539   f = test_LU_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
 540   gsl_test(f, "  LU_solve vander(3)");
 541   s += f;
 542
 543   f = test_LU_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
 544   gsl_test(f, "  LU_solve vander(4)");
 545   s += f;
 546
 547   f = test_LU_solve_dim(vander12, vander12_solution, 0.05);
 548   gsl_test(f, "  LU_solve vander(12)");
 549   s += f;
 550
 551   return s;
 552 }
 553
 554
 555 int
 556 test_LUc_solve_dim(const gsl_matrix_complex * m, const double * actual, double eps)
 557 {
 558   int s = 0;
 559   int signum;
 560   unsigned long i, dim = m->size1;
 561
 562   gsl_permutation * perm = gsl_permutation_alloc(dim);
 563   gsl_vector_complex * rhs = gsl_vector_complex_alloc(dim);
 564   gsl_matrix_complex * lu  = gsl_matrix_complex_alloc(dim,dim);
 565   gsl_vector_complex * x = gsl_vector_complex_alloc(dim);
 566   gsl_vector_complex * residual = gsl_vector_complex_alloc(dim);
 567   gsl_matrix_complex_memcpy(lu,m);
 568   for(i=0; i<dim; i++)
 569     {
 570       gsl_complex z = gsl_complex_rect (2.0*i+1.0, 2.0*i+2.0);
 571       gsl_vector_complex_set(rhs, i, z);
 572     }
 573   s += gsl_linalg_complex_LU_decomp(lu, perm, &signum);
 574   s += gsl_linalg_complex_LU_solve(lu, perm, rhs, x);
 575
 576   for(i=0; i<dim; i++) {
 577     gsl_complex z = gsl_vector_complex_get(x, i);
 578     int foo_r = check(GSL_REAL(z),actual[2*i],eps);
 579     int foo_i = check(GSL_IMAG(z),actual[2*i+1],eps);
 580     if(foo_r || foo_i) {
 581       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, GSL_REAL(z), actual[2*i]);
 582       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, GSL_IMAG(z), actual[2*i+1]);
 583     }
 584     s += foo_r + foo_i;
 585   }
 586
 587   s += gsl_linalg_complex_LU_refine(m, lu, perm, rhs, x, residual);
 588
 589   for(i=0; i<dim; i++) {
 590     gsl_complex z = gsl_vector_complex_get(x, i);
 591     int foo_r = check(GSL_REAL(z),actual[2*i],eps);
 592     int foo_i = check(GSL_IMAG(z),actual[2*i+1],eps);
 593     if(foo_r || foo_i) {
 594       printf("%3lu[%lu]: %22.18g   %22.18g (improved)\n", dim, i, GSL_REAL(z), actual[2*i]);
 595       printf("%3lu[%lu]: %22.18g   %22.18g (improved)\n", dim, i, GSL_IMAG(z), actual[2*i+1]);
 596     }
 597     s += foo_r + foo_i;
 598   }
 599
 600   gsl_vector_complex_free(residual);
 601   gsl_vector_complex_free(x);
 602   gsl_matrix_complex_free(lu);
 603   gsl_vector_complex_free(rhs);
 604   gsl_permutation_free(perm);
 605
 606   return s;
 607 }
 608
 609
 610 int test_LUc_solve(void)
 611 {
 612   int f;
 613   int s = 0;
 614
 615   f = test_LUc_solve_dim(c7, c7_solution, 1024.0 * 1024.0 * GSL_DBL_EPSILON);
 616   gsl_test(f, "  complex_LU_solve complex(7)");
 617   s += f;
 618
 619   return s;
 620 }
 621
 622
 623 int
 624 test_QR_solve_dim(const gsl_matrix * m, const double * actual, double eps)
 625 {
 626   int s = 0;
 627   unsigned long i, dim = m->size1;
 628
 629   gsl_vector * rhs = gsl_vector_alloc(dim);
 630   gsl_matrix * qr  = gsl_matrix_alloc(dim,dim);
 631   gsl_vector * d = gsl_vector_alloc(dim);
 632   gsl_vector * x = gsl_vector_alloc(dim);
 633
 634   gsl_matrix_memcpy(qr,m);
 635   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
 636   s += gsl_linalg_QR_decomp(qr, d);
 637   s += gsl_linalg_QR_solve(qr, d, rhs, x);
 638   for(i=0; i<dim; i++) {
 639     int foo = check(gsl_vector_get(x, i), actual[i], eps);
 640     if(foo) {
 641       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
 642     }
 643     s += foo;
 644   }
 645
 646   gsl_vector_free(x);
 647   gsl_vector_free(d);
 648   gsl_matrix_free(qr);
 649   gsl_vector_free(rhs);
 650
 651   return s;
 652 }
 653
 654 int test_QR_solve(void)
 655 {
 656   int f;
 657   int s = 0;
 658
 659   f = test_QR_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
 660   gsl_test(f, "  QR_solve hilbert(2)");
 661   s += f;
 662
 663   f = test_QR_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
 664   gsl_test(f, "  QR_solve hilbert(3)");
 665   s += f;
 666
 667   f = test_QR_solve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
 668   gsl_test(f, "  QR_solve hilbert(4)");
 669   s += f;
 670
 671   f = test_QR_solve_dim(hilb12, hilb12_solution, 0.5);
 672   gsl_test(f, "  QR_solve hilbert(12)");
 673   s += f;
 674
 675   f = test_QR_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
 676   gsl_test(f, "  QR_solve vander(2)");
 677   s += f;
 678
 679   f = test_QR_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
 680   gsl_test(f, "  QR_solve vander(3)");
 681   s += f;
 682
 683   f = test_QR_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
 684   gsl_test(f, "  QR_solve vander(4)");
 685   s += f;
 686
 687   f = test_QR_solve_dim(vander12, vander12_solution, 0.05);
 688   gsl_test(f, "  QR_solve vander(12)");
 689   s += f;
 690
 691   return s;
 692 }
 693
 694
 695 int
 696 test_QR_QRsolve_dim(const gsl_matrix * m, const double * actual, double eps)
 697 {
 698   int s = 0;
 699   unsigned long i, dim = m->size1;
 700
 701   gsl_vector * rhs = gsl_vector_alloc(dim);
 702   gsl_matrix * qr  = gsl_matrix_alloc(dim,dim);
 703   gsl_matrix * q  = gsl_matrix_alloc(dim,dim);
 704   gsl_matrix * r  = gsl_matrix_alloc(dim,dim);
 705   gsl_vector * d = gsl_vector_alloc(dim);
 706   gsl_vector * x = gsl_vector_alloc(dim);
 707
 708   gsl_matrix_memcpy(qr,m);
 709   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
 710   s += gsl_linalg_QR_decomp(qr, d);
 711   s += gsl_linalg_QR_unpack(qr, d, q, r);
 712   s += gsl_linalg_QR_QRsolve(q, r, rhs, x);
 713   for(i=0; i<dim; i++) {
 714     int foo = check(gsl_vector_get(x, i), actual[i], eps);
 715     if(foo) {
 716       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
 717     }
 718     s += foo;
 719   }
 720
 721   gsl_vector_free(x);
 722   gsl_vector_free(d);
 723   gsl_matrix_free(qr);
 724   gsl_matrix_free(q);
 725   gsl_matrix_free(r);
 726   gsl_vector_free(rhs);
 727   return s;
 728 }
 729
 730 int test_QR_QRsolve(void)
 731 {
 732   int f;
 733   int s = 0;
 734
 735   f = test_QR_QRsolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
 736   gsl_test(f, "  QR_QRsolve hilbert(2)");
 737   s += f;
 738
 739   f = test_QR_QRsolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
 740   gsl_test(f, "  QR_QRsolve hilbert(3)");
 741   s += f;
 742
 743   f = test_QR_QRsolve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
 744   gsl_test(f, "  QR_QRsolve hilbert(4)");
 745   s += f;
 746
 747   f = test_QR_QRsolve_dim(hilb12, hilb12_solution, 0.5);
 748   gsl_test(f, "  QR_QRsolve hilbert(12)");
 749   s += f;
 750
 751   f = test_QR_QRsolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
 752   gsl_test(f, "  QR_QRsolve vander(2)");
 753   s += f;
 754
 755   f = test_QR_QRsolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
 756   gsl_test(f, "  QR_QRsolve vander(3)");
 757   s += f;
 758
 759   f = test_QR_QRsolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
 760   gsl_test(f, "  QR_QRsolve vander(4)");
 761   s += f;
 762
 763   f = test_QR_QRsolve_dim(vander12, vander12_solution, 0.05);
 764   gsl_test(f, "  QR_QRsolve vander(12)");
 765   s += f;
 766
 767   return s;
 768 }
 769
 770
 771 int
 772 test_QR_lssolve_dim(const gsl_matrix * m, const double * actual, double eps)
 773 {
 774   int s = 0;
 775   unsigned long i, M = m->size1, N = m->size2;
 776
 777   gsl_vector * rhs = gsl_vector_alloc(M);
 778   gsl_matrix * qr  = gsl_matrix_alloc(M,N);
 779   gsl_vector * d = gsl_vector_alloc(N);
 780   gsl_vector * x = gsl_vector_alloc(N);
 781   gsl_vector * r = gsl_vector_alloc(M);
 782   gsl_vector * res = gsl_vector_alloc(M);
 783
 784   gsl_matrix_memcpy(qr,m);
 785   for(i=0; i<M; i++) gsl_vector_set(rhs, i, i+1.0);
 786   s += gsl_linalg_QR_decomp(qr, d);
 787   s += gsl_linalg_QR_lssolve(qr, d, rhs, x, res);
 788
 789   for(i=0; i<N; i++) {
 790     int foo = check(gsl_vector_get(x, i), actual[i], eps);
 791     if(foo) {
 792       printf("(%3lu,%3lu)[%lu]: %22.18g   %22.18g\n", M, N, i, gsl_vector_get(x, i), actual[i]);
 793     }
 794     s += foo;
 795   }
 796
 797   /* compute residual r = b - m x */
 798   if (M == N) {
 799     gsl_vector_set_zero(r);
 800   } else {
 801     gsl_vector_memcpy(r, rhs);
 802     gsl_blas_dgemv(CblasNoTrans, -1.0, m, x, 1.0, r);
 803   };
 804
 805   for(i=0; i<N; i++) {
 806     int foo = check(gsl_vector_get(res, i), gsl_vector_get(r,i), sqrt(eps));
 807     if(foo) {
 808       printf("(%3lu,%3lu)[%lu]: %22.18g   %22.18g\n", M, N, i, gsl_vector_get(res, i), gsl_vector_get(r,i));
 809     }
 810     s += foo;
 811   }
 812
 813   gsl_vector_free(r);
 814   gsl_vector_free(res);
 815   gsl_vector_free(x);
 816   gsl_vector_free(d);
 817   gsl_matrix_free(qr);
 818   gsl_vector_free(rhs);
 819
 820   return s;
 821 }
 822
 823 int test_QR_lssolve(void)
 824 {
 825   int f;
 826   int s = 0;
 827
 828   f = test_QR_lssolve_dim(m53, m53_lssolution, 2 * 64.0 * GSL_DBL_EPSILON);
 829   gsl_test(f, "  QR_lssolve m(5,3)");
 830   s += f;
 831
 832   f = test_QR_lssolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
 833   gsl_test(f, "  QR_lssolve hilbert(2)");
 834   s += f;
 835
 836   f = test_QR_lssolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
 837   gsl_test(f, "  QR_lssolve hilbert(3)");
 838   s += f;
 839
 840   f = test_QR_lssolve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
 841   gsl_test(f, "  QR_lssolve hilbert(4)");
 842   s += f;
 843
 844   f = test_QR_lssolve_dim(hilb12, hilb12_solution, 0.5);
 845   gsl_test(f, "  QR_lssolve hilbert(12)");
 846   s += f;
 847
 848   f = test_QR_lssolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
 849   gsl_test(f, "  QR_lssolve vander(2)");
 850   s += f;
 851
 852   f = test_QR_lssolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
 853   gsl_test(f, "  QR_lssolve vander(3)");
 854   s += f;
 855
 856   f = test_QR_lssolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
 857   gsl_test(f, "  QR_lssolve vander(4)");
 858   s += f;
 859
 860   f = test_QR_lssolve_dim(vander12, vander12_solution, 0.05);
 861   gsl_test(f, "  QR_lssolve vander(12)");
 862   s += f;
 863
 864   return s;
 865 }
 866
 867
 868 int
 869 test_QR_decomp_dim(const gsl_matrix * m, double eps)
 870 {
 871   int s = 0;
 872   unsigned long i,j, M = m->size1, N = m->size2;
 873
 874   gsl_matrix * qr = gsl_matrix_alloc(M,N);
 875   gsl_matrix * a  = gsl_matrix_alloc(M,N);
 876   gsl_matrix * q  = gsl_matrix_alloc(M,M);
 877   gsl_matrix * r  = gsl_matrix_alloc(M,N);
 878   gsl_vector * d = gsl_vector_alloc(GSL_MIN(M,N));
 879
 880   gsl_matrix_memcpy(qr,m);
 881
 882   s += gsl_linalg_QR_decomp(qr, d);
 883   s += gsl_linalg_QR_unpack(qr, d, q, r);
 884
 885   /* compute a = q r */
 886   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, q, r, 0.0, a);
 887
 888   for(i=0; i<M; i++) {
 889     for(j=0; j<N; j++) {
 890       double aij = gsl_matrix_get(a, i, j);
 891       double mij = gsl_matrix_get(m, i, j);
 892       int foo = check(aij, mij, eps);
 893       if(foo) {
 894         printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
 895       }
 896       s += foo;
 897     }
 898   }
 899
 900   gsl_vector_free(d);
 901   gsl_matrix_free(qr);
 902   gsl_matrix_free(a);
 903   gsl_matrix_free(q);
 904   gsl_matrix_free(r);
 905
 906   return s;
 907 }
 908
 909 int test_QR_decomp(void)
 910 {
 911   int f;
 912   int s = 0;
 913
 914   f = test_QR_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
 915   gsl_test(f, "  QR_decomp m(3,5)");
 916   s += f;
 917
 918   f = test_QR_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
 919   gsl_test(f, "  QR_decomp m(5,3)");
 920   s += f;
 921
 922   f = test_QR_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
 923   gsl_test(f, "  QR_decomp hilbert(2)");
 924   s += f;
 925
 926   f = test_QR_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
 927   gsl_test(f, "  QR_decomp hilbert(3)");
 928   s += f;
 929
 930   f = test_QR_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
 931   gsl_test(f, "  QR_decomp hilbert(4)");
 932   s += f;
 933
 934   f = test_QR_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
 935   gsl_test(f, "  QR_decomp hilbert(12)");
 936   s += f;
 937
 938   f = test_QR_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
 939   gsl_test(f, "  QR_decomp vander(2)");
 940   s += f;
 941
 942   f = test_QR_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
 943   gsl_test(f, "  QR_decomp vander(3)");
 944   s += f;
 945
 946   f = test_QR_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
 947   gsl_test(f, "  QR_decomp vander(4)");
 948   s += f;
 949
 950   f = test_QR_decomp_dim(vander12, 0.0005); /* FIXME: bad accuracy */
 951   gsl_test(f, "  QR_decomp vander(12)");
 952   s += f;
 953
 954   return s;
 955 }
 956
 957 int
 958 test_QRPT_solve_dim(const gsl_matrix * m, const double * actual, double eps)
 959 {
 960   int s = 0;
 961   int signum;
 962   unsigned long i, dim = m->size1;
 963
 964   gsl_permutation * perm = gsl_permutation_alloc(dim);
 965   gsl_vector * rhs = gsl_vector_alloc(dim);
 966   gsl_matrix * qr  = gsl_matrix_alloc(dim,dim);
 967   gsl_vector * d = gsl_vector_alloc(dim);
 968   gsl_vector * x = gsl_vector_alloc(dim);
 969   gsl_vector * norm = gsl_vector_alloc(dim);
 970
 971   gsl_matrix_memcpy(qr,m);
 972   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
 973   s += gsl_linalg_QRPT_decomp(qr, d, perm, &signum, norm);
 974   s += gsl_linalg_QRPT_solve(qr, d, perm, rhs, x);
 975   for(i=0; i<dim; i++) {
 976     int foo = check(gsl_vector_get(x, i), actual[i], eps);
 977     if(foo) {
 978       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
 979     }
 980     s += foo;
 981   }
 982
 983   gsl_vector_free(norm);
 984   gsl_vector_free(x);
 985   gsl_vector_free(d);
 986   gsl_matrix_free(qr);
 987   gsl_vector_free(rhs);
 988   gsl_permutation_free(perm);
 989
 990   return s;
 991 }
 992
 993 int test_QRPT_solve(void)
 994 {
 995   int f;
 996   int s = 0;
 997
 998   f = test_QRPT_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
 999   gsl_test(f, "  QRPT_solve hilbert(2)");
1000   s += f;
1001
1002   f = test_QRPT_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
1003   gsl_test(f, "  QRPT_solve hilbert(3)");
1004   s += f;
1005
1006   f = test_QRPT_solve_dim(hilb4, hilb4_solution, 2 * 2048.0 * GSL_DBL_EPSILON);
1007   gsl_test(f, "  QRPT_solve hilbert(4)");
1008   s += f;
1009
1010   f = test_QRPT_solve_dim(hilb12, hilb12_solution, 0.5);
1011   gsl_test(f, "  QRPT_solve hilbert(12)");
1012   s += f;
1013
1014   f = test_QRPT_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
1015   gsl_test(f, "  QRPT_solve vander(2)");
1016   s += f;
1017
1018   f = test_QRPT_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
1019   gsl_test(f, "  QRPT_solve vander(3)");
1020   s += f;
1021
1022   f = test_QRPT_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
1023   gsl_test(f, "  QRPT_solve vander(4)");
1024   s += f;
1025
1026   f = test_QRPT_solve_dim(vander12, vander12_solution, 0.05);
1027   gsl_test(f, "  QRPT_solve vander(12)");
1028   s += f;
1029
1030   return s;
1031 }
1032
1033 int
1034 test_QRPT_QRsolve_dim(const gsl_matrix * m, const double * actual, double eps)
1035 {
1036   int s = 0;
1037   int signum;
1038   unsigned long i, dim = m->size1;
1039
1040   gsl_permutation * perm = gsl_permutation_alloc(dim);
1041   gsl_vector * rhs = gsl_vector_alloc(dim);
1042   gsl_matrix * qr  = gsl_matrix_alloc(dim,dim);
1043   gsl_matrix * q  = gsl_matrix_alloc(dim,dim);
1044   gsl_matrix * r  = gsl_matrix_alloc(dim,dim);
1045   gsl_vector * d = gsl_vector_alloc(dim);
1046   gsl_vector * x = gsl_vector_alloc(dim);
1047   gsl_vector * norm = gsl_vector_alloc(dim);
1048
1049   gsl_matrix_memcpy(qr,m);
1050   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
1051   s += gsl_linalg_QRPT_decomp2(qr, q, r, d, perm, &signum, norm);
1052   s += gsl_linalg_QRPT_QRsolve(q, r, perm, rhs, x);
1053   for(i=0; i<dim; i++) {
1054     int foo = check(gsl_vector_get(x, i), actual[i], eps);
1055     if(foo) {
1056       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
1057     }
1058     s += foo;
1059   }
1060
1061   gsl_vector_free(norm);
1062   gsl_vector_free(x);
1063   gsl_vector_free(d);
1064   gsl_matrix_free(qr);
1065   gsl_matrix_free(q);
1066   gsl_matrix_free(r);
1067   gsl_vector_free(rhs);
1068   gsl_permutation_free(perm);
1069
1070   return s;
1071 }
1072
1073 int test_QRPT_QRsolve(void)
1074 {
1075   int f;
1076   int s = 0;
1077
1078   f = test_QRPT_QRsolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
1079   gsl_test(f, "  QRPT_QRsolve hilbert(2)");
1080   s += f;
1081
1082   f = test_QRPT_QRsolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
1083   gsl_test(f, "  QRPT_QRsolve hilbert(3)");
1084   s += f;
1085
1086   f = test_QRPT_QRsolve_dim(hilb4, hilb4_solution, 2 * 2048.0 * GSL_DBL_EPSILON);
1087   gsl_test(f, "  QRPT_QRsolve hilbert(4)");
1088   s += f;
1089
1090   f = test_QRPT_QRsolve_dim(hilb12, hilb12_solution, 0.5);
1091   gsl_test(f, "  QRPT_QRsolve hilbert(12)");
1092   s += f;
1093
1094   f = test_QRPT_QRsolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
1095   gsl_test(f, "  QRPT_QRsolve vander(2)");
1096   s += f;
1097
1098   f = test_QRPT_QRsolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
1099   gsl_test(f, "  QRPT_QRsolve vander(3)");
1100   s += f;
1101
1102   f = test_QRPT_QRsolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
1103   gsl_test(f, "  QRPT_QRsolve vander(4)");
1104   s += f;
1105
1106   f = test_QRPT_QRsolve_dim(vander12, vander12_solution, 0.05);
1107   gsl_test(f, "  QRPT_QRsolve vander(12)");
1108   s += f;
1109
1110   return s;
1111 }
1112
1113 int
1114 test_QRPT_decomp_dim(const gsl_matrix * m, double eps)
1115 {
1116   int s = 0, signum;
1117   unsigned long i,j, M = m->size1, N = m->size2;
1118
1119   gsl_matrix * qr = gsl_matrix_alloc(M,N);
1120   gsl_matrix * a  = gsl_matrix_alloc(M,N);
1121   gsl_matrix * q  = gsl_matrix_alloc(M,M);
1122   gsl_matrix * r  = gsl_matrix_alloc(M,N);
1123   gsl_vector * d = gsl_vector_alloc(GSL_MIN(M,N));
1124   gsl_vector * norm = gsl_vector_alloc(N);
1125
1126   gsl_permutation * perm = gsl_permutation_alloc(N);
1127
1128   gsl_matrix_memcpy(qr,m);
1129
1130   s += gsl_linalg_QRPT_decomp(qr, d, perm, &signum, norm);
1131   s += gsl_linalg_QR_unpack(qr, d, q, r);
1132
1133   /* compute a = q r */
1134   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, q, r, 0.0, a);
1135
1136
1137   /* Compute QR P^T by permuting the elements of the rows of QR */
1138
1139   for (i = 0; i < M; i++) {
1140     gsl_vector_view row = gsl_matrix_row (a, i);
1141     gsl_permute_vector_inverse (perm, &row.vector);
1142   }
1143
1144   for(i=0; i<M; i++) {
1145     for(j=0; j<N; j++) {
1146       double aij = gsl_matrix_get(a, i, j);
1147       double mij = gsl_matrix_get(m, i, j);
1148       int foo = check(aij, mij, eps);
1149       if(foo) {
1150         printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
1151       }
1152       s += foo;
1153     }
1154   }
1155
1156   gsl_permutation_free (perm);
1157   gsl_vector_free(norm);
1158   gsl_vector_free(d);
1159   gsl_matrix_free(qr);
1160   gsl_matrix_free(a);
1161   gsl_matrix_free(q);
1162   gsl_matrix_free(r);
1163
1164   return s;
1165 }
1166
1167 int test_QRPT_decomp(void)
1168 {
1169   int f;
1170   int s = 0;
1171
1172   f = test_QRPT_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
1173   gsl_test(f, "  QRPT_decomp m(3,5)");
1174   s += f;
1175
1176   f = test_QRPT_decomp_dim(m53, 2 * 8.0 * GSL_DBL_EPSILON);
1177   gsl_test(f, "  QRPT_decomp m(5,3)");
1178   s += f;
1179
1180   f = test_QRPT_decomp_dim(s35, 2 * 8.0 * GSL_DBL_EPSILON);
1181   gsl_test(f, "  QRPT_decomp s(3,5)");
1182   s += f;
1183
1184   f = test_QRPT_decomp_dim(s53, 2 * 8.0 * GSL_DBL_EPSILON);
1185   gsl_test(f, "  QRPT_decomp s(5,3)");
1186   s += f;
1187
1188   f = test_QRPT_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
1189   gsl_test(f, "  QRPT_decomp hilbert(2)");
1190   s += f;
1191
1192   f = test_QRPT_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
1193   gsl_test(f, "  QRPT_decomp hilbert(3)");
1194   s += f;
1195
1196   f = test_QRPT_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
1197   gsl_test(f, "  QRPT_decomp hilbert(4)");
1198   s += f;
1199
1200   f = test_QRPT_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
1201   gsl_test(f, "  QRPT_decomp hilbert(12)");
1202   s += f;
1203
1204   f = test_QRPT_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
1205   gsl_test(f, "  QRPT_decomp vander(2)");
1206   s += f;
1207
1208   f = test_QRPT_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
1209   gsl_test(f, "  QRPT_decomp vander(3)");
1210   s += f;
1211
1212   f = test_QRPT_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
1213   gsl_test(f, "  QRPT_decomp vander(4)");
1214   s += f;
1215
1216   f = test_QRPT_decomp_dim(vander12, 0.0005); /* FIXME: bad accuracy */
1217   gsl_test(f, "  QRPT_decomp vander(12)");
1218   s += f;
1219
1220   return s;
1221 }
1222
1223
1224 int
1225 test_QR_update_dim(const gsl_matrix * m, double eps)
1226 {
1227   int s = 0;
1228   unsigned long i,j,k, M = m->size1, N = m->size2;
1229
1230   gsl_vector * rhs = gsl_vector_alloc(N);
1231   gsl_matrix * qr1  = gsl_matrix_alloc(M,N);
1232   gsl_matrix * qr2  = gsl_matrix_alloc(M,N);
1233   gsl_matrix * q1  = gsl_matrix_alloc(M,M);
1234   gsl_matrix * r1  = gsl_matrix_alloc(M,N);
1235   gsl_matrix * q2  = gsl_matrix_alloc(M,M);
1236   gsl_matrix * r2  = gsl_matrix_alloc(M,N);
1237   gsl_vector * d = gsl_vector_alloc(GSL_MIN(M,N));
1238   gsl_vector * solution1 = gsl_vector_alloc(N);
1239   gsl_vector * solution2 = gsl_vector_alloc(N);
1240   gsl_vector * u = gsl_vector_alloc(M);
1241   gsl_vector * v = gsl_vector_alloc(N);
1242   gsl_vector * w = gsl_vector_alloc(M);
1243
1244   gsl_matrix_memcpy(qr1,m);
1245   gsl_matrix_memcpy(qr2,m);
1246   for(i=0; i<N; i++) gsl_vector_set(rhs, i, i+1.0);
1247   for(i=0; i<M; i++) gsl_vector_set(u, i, sin(i+1.0));
1248   for(i=0; i<N; i++) gsl_vector_set(v, i, cos(i+2.0) + sin(i*i+3.0));
1249
1250   for(i=0; i<M; i++)
1251     {
1252       double ui = gsl_vector_get(u, i);
1253       for(j=0; j<N; j++)
1254         {
1255           double vj = gsl_vector_get(v, j);
1256           double qij = gsl_matrix_get(qr1, i, j);
1257           gsl_matrix_set(qr1, i, j, qij + ui * vj);
1258         }
1259     }
1260
1261   s += gsl_linalg_QR_decomp(qr2, d);
1262   s += gsl_linalg_QR_unpack(qr2, d, q2, r2);
1263
1264   /* compute w = Q^T u */
1265
1266   for (j = 0; j < M; j++)
1267     {
1268       double sum = 0;
1269       for (i = 0; i < M; i++)
1270           sum += gsl_matrix_get (q2, i, j) * gsl_vector_get (u, i);
1271       gsl_vector_set (w, j, sum);
1272     }
1273
1274   s += gsl_linalg_QR_update(q2, r2, w, v);
1275
1276   /* compute qr2 = q2 * r2 */
1277
1278   for (i = 0; i < M; i++)
1279     {
1280       for (j = 0; j< N; j++)
1281         {
1282           double sum = 0;
1283           for (k = 0; k <= GSL_MIN(j,M-1); k++)
1284             {
1285               double qik = gsl_matrix_get(q2, i, k);
1286               double rkj = gsl_matrix_get(r2, k, j);
1287               sum += qik * rkj ;
1288             }
1289           gsl_matrix_set (qr2, i, j, sum);
1290         }
1291     }
1292
1293   for(i=0; i<M; i++) {
1294     for(j=0; j<N; j++) {
1295       double s1 = gsl_matrix_get(qr1, i, j);
1296       double s2 = gsl_matrix_get(qr2, i, j);
1297
1298       int foo = check(s1, s2, eps);
1299       if(foo) {
1300         printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, s1, s2);
1301       }
1302       s += foo;
1303     }
1304   }
1305
1306   gsl_vector_free(solution1);
1307   gsl_vector_free(solution2);
1308   gsl_vector_free(d);
1309   gsl_vector_free(u);
1310   gsl_vector_free(v);
1311   gsl_vector_free(w);
1312   gsl_matrix_free(qr1);
1313   gsl_matrix_free(qr2);
1314   gsl_matrix_free(q1);
1315   gsl_matrix_free(r1);
1316   gsl_matrix_free(q2);
1317   gsl_matrix_free(r2);
1318   gsl_vector_free(rhs);
1319
1320   return s;
1321 }
1322
1323 int test_QR_update(void)
1324 {
1325   int f;
1326   int s = 0;
1327
1328   f = test_QR_update_dim(m35, 2 * 512.0 * GSL_DBL_EPSILON);
1329   gsl_test(f, "  QR_update m(3,5)");
1330   s += f;
1331
1332   f = test_QR_update_dim(m53, 2 * 512.0 * GSL_DBL_EPSILON);
1333   gsl_test(f, "  QR_update m(5,3)");
1334   s += f;
1335
1336   f = test_QR_update_dim(hilb2,  2 * 512.0 * GSL_DBL_EPSILON);
1337   gsl_test(f, "  QR_update hilbert(2)");
1338   s += f;
1339
1340   f = test_QR_update_dim(hilb3,  2 * 512.0 * GSL_DBL_EPSILON);
1341   gsl_test(f, "  QR_update hilbert(3)");
1342   s += f;
1343
1344   f = test_QR_update_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
1345   gsl_test(f, "  QR_update hilbert(4)");
1346   s += f;
1347
1348   f = test_QR_update_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
1349   gsl_test(f, "  QR_update hilbert(12)");
1350   s += f;
1351
1352   f = test_QR_update_dim(vander2, 8.0 * GSL_DBL_EPSILON);
1353   gsl_test(f, "  QR_update vander(2)");
1354   s += f;
1355
1356   f = test_QR_update_dim(vander3, 64.0 * GSL_DBL_EPSILON);
1357   gsl_test(f, "  QR_update vander(3)");
1358   s += f;
1359
1360   f = test_QR_update_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
1361   gsl_test(f, "  QR_update vander(4)");
1362   s += f;
1363
1364   f = test_QR_update_dim(vander12, 0.0005); /* FIXME: bad accuracy */
1365   gsl_test(f, "  QR_update vander(12)");
1366   s += f;
1367
1368   return s;
1369 }
1370
1371 int
1372 test_LQ_solve_dim(const gsl_matrix * m, const double * actual, double eps)
1373 {
1374   int s = 0;
1375   unsigned long i, dim = m->size1;
1376
1377   gsl_vector * rhs = gsl_vector_alloc(dim);
1378   gsl_matrix * lq  = gsl_matrix_alloc(dim,dim);
1379   gsl_vector * d = gsl_vector_alloc(dim);
1380   gsl_vector * x = gsl_vector_alloc(dim);
1381
1382   gsl_matrix_transpose_memcpy(lq,m);
1383   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
1384   s += gsl_linalg_LQ_decomp(lq, d);
1385   s += gsl_linalg_LQ_solve_T(lq, d, rhs, x);
1386   for(i=0; i<dim; i++) {
1387     int foo = check(gsl_vector_get(x, i), actual[i], eps);
1388     if(foo) {
1389       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
1390     }
1391     s += foo;
1392   }
1393
1394   gsl_vector_free(x);
1395   gsl_vector_free(d);
1396   gsl_matrix_free(lq);
1397   gsl_vector_free(rhs);
1398
1399   return s;
1400 }
1401
1402 int test_LQ_solve(void)
1403 {
1404   int f;
1405   int s = 0;
1406
1407   f = test_LQ_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
1408   gsl_test(f, "  LQ_solve hilbert(2)");
1409   s += f;
1410
1411   f = test_LQ_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
1412   gsl_test(f, "  LQ_solve hilbert(3)");
1413   s += f;
1414
1415   f = test_LQ_solve_dim(hilb4, hilb4_solution, 4 * 1024.0 * GSL_DBL_EPSILON);
1416   gsl_test(f, "  LQ_solve hilbert(4)");
1417   s += f;
1418
1419   f = test_LQ_solve_dim(hilb12, hilb12_solution, 0.5);
1420   gsl_test(f, "  LQ_solve hilbert(12)");
1421   s += f;
1422
1423   f = test_LQ_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
1424   gsl_test(f, "  LQ_solve vander(2)");
1425   s += f;
1426
1427   f = test_LQ_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
1428   gsl_test(f, "  LQ_solve vander(3)");
1429   s += f;
1430
1431   f = test_LQ_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
1432   gsl_test(f, "  LQ_solve vander(4)");
1433   s += f;
1434
1435   f = test_LQ_solve_dim(vander12, vander12_solution, 0.05);
1436   gsl_test(f, "  LQ_solve vander(12)");
1437   s += f;
1438
1439   return s;
1440 }
1441
1442
1443
1444
1445 int
1446 test_LQ_LQsolve_dim(const gsl_matrix * m, const double * actual, double eps)
1447 {
1448   int s = 0;
1449   unsigned long i, dim = m->size1;
1450
1451   gsl_vector * rhs = gsl_vector_alloc(dim);
1452   gsl_matrix * lq  = gsl_matrix_alloc(dim,dim);
1453   gsl_matrix * q  = gsl_matrix_alloc(dim,dim);
1454   gsl_matrix * l  = gsl_matrix_alloc(dim,dim);
1455   gsl_vector * d = gsl_vector_alloc(dim);
1456   gsl_vector * x = gsl_vector_alloc(dim);
1457
1458   gsl_matrix_transpose_memcpy(lq,m);
1459   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
1460   s += gsl_linalg_LQ_decomp(lq, d);
1461   s += gsl_linalg_LQ_unpack(lq, d, q, l);
1462   s += gsl_linalg_LQ_LQsolve(q, l, rhs, x);
1463   for(i=0; i<dim; i++) {
1464     int foo = check(gsl_vector_get(x, i), actual[i], eps);
1465     if(foo) {
1466       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
1467     }
1468     s += foo;
1469   }
1470
1471   gsl_vector_free(x);
1472   gsl_vector_free(d);
1473   gsl_matrix_free(lq);
1474   gsl_matrix_free(q);
1475   gsl_matrix_free(l);
1476   gsl_vector_free(rhs);
1477
1478   return s;
1479 }
1480
1481 int test_LQ_LQsolve(void)
1482 {
1483   int f;
1484   int s = 0;
1485
1486   f = test_LQ_LQsolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
1487   gsl_test(f, "  LQ_LQsolve hilbert(2)");
1488   s += f;
1489
1490   f = test_LQ_LQsolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
1491   gsl_test(f, "  LQ_LQsolve hilbert(3)");
1492   s += f;
1493
1494   f = test_LQ_LQsolve_dim(hilb4, hilb4_solution, 4 * 1024.0 * GSL_DBL_EPSILON);
1495   gsl_test(f, "  LQ_LQsolve hilbert(4)");
1496   s += f;
1497
1498   f = test_LQ_LQsolve_dim(hilb12, hilb12_solution, 0.5);
1499   gsl_test(f, "  LQ_LQsolve hilbert(12)");
1500   s += f;
1501
1502   f = test_LQ_LQsolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
1503   gsl_test(f, "  LQ_LQsolve vander(2)");
1504   s += f;
1505
1506   f = test_LQ_LQsolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
1507   gsl_test(f, "  LQ_LQsolve vander(3)");
1508   s += f;
1509
1510   f = test_LQ_LQsolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
1511   gsl_test(f, "  LQ_LQsolve vander(4)");
1512   s += f;
1513
1514   f = test_LQ_LQsolve_dim(vander12, vander12_solution, 0.05);
1515   gsl_test(f, "  LQ_LQsolve vander(12)");
1516   s += f;
1517
1518   return s;
1519 }
1520
1521
1522 int
1523 test_LQ_lssolve_dim(const gsl_matrix * m, const double * actual, double eps)
1524 {
1525   int s = 0;
1526   unsigned long i, M = m->size1, N = m->size2;
1527
1528   gsl_vector * rhs = gsl_vector_alloc(M);
1529   gsl_matrix * lq  = gsl_matrix_alloc(N,M);
1530   gsl_vector * d = gsl_vector_alloc(N);
1531   gsl_vector * x = gsl_vector_alloc(N);
1532   gsl_vector * r = gsl_vector_alloc(M);
1533   gsl_vector * res = gsl_vector_alloc(M);
1534
1535   gsl_matrix_transpose_memcpy(lq,m);
1536   for(i=0; i<M; i++) gsl_vector_set(rhs, i, i+1.0);
1537   s += gsl_linalg_LQ_decomp(lq, d);
1538   s += gsl_linalg_LQ_lssolve_T(lq, d, rhs, x, res);
1539
1540   for(i=0; i<N; i++) {
1541     int foo = check(gsl_vector_get(x, i), actual[i], eps);
1542     if(foo) {
1543       printf("(%3lu,%3lu)[%lu]: %22.18g   %22.18g\n", M, N, i, gsl_vector_get(x, i), actual[i]);
1544     }
1545     s += foo;
1546   }
1547
1548
1549    /* compute residual r = b - m x */
1550   if (M == N) {
1551     gsl_vector_set_zero(r);
1552   } else {
1553     gsl_vector_memcpy(r, rhs);
1554     gsl_blas_dgemv(CblasNoTrans, -1.0, m, x, 1.0, r);
1555   };
1556
1557   for(i=0; i<N; i++) {
1558     int foo = check(gsl_vector_get(res, i), gsl_vector_get(r,i), sqrt(eps));
1559     if(foo) {
1560       printf("(%3lu,%3lu)[%lu]: %22.18g   %22.18g\n", M, N, i, gsl_vector_get(res, i), gsl_vector_get(r,i));
1561     }
1562     s += foo;
1563   }
1564
1565   gsl_vector_free(r);
1566   gsl_vector_free(res);
1567   gsl_vector_free(x);
1568   gsl_vector_free(d);
1569   gsl_matrix_free(lq);
1570   gsl_vector_free(rhs);
1571
1572   return s;
1573 }
1574
1575 int test_LQ_lssolve(void)
1576 {
1577   int f;
1578   int s = 0;
1579
1580   f = test_LQ_lssolve_dim(m53, m53_lssolution, 2 * 64.0 * GSL_DBL_EPSILON);
1581   gsl_test(f, "  LQ_lssolve m(5,3)");
1582   s += f;
1583
1584   f = test_LQ_lssolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
1585   gsl_test(f, "  LQ_lssolve hilbert(2)");
1586   s += f;
1587
1588   f = test_LQ_lssolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
1589   gsl_test(f, "  LQ_lssolve hilbert(3)");
1590   s += f;
1591
1592   f = test_LQ_lssolve_dim(hilb4, hilb4_solution, 4 * 1024.0 * GSL_DBL_EPSILON);
1593   gsl_test(f, "  LQ_lssolve hilbert(4)");
1594   s += f;
1595
1596   f = test_LQ_lssolve_dim(hilb12, hilb12_solution, 0.5);
1597   gsl_test(f, "  LQ_lssolve hilbert(12)");
1598   s += f;
1599
1600   f = test_LQ_lssolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
1601   gsl_test(f, "  LQ_lssolve vander(2)");
1602   s += f;
1603
1604   f = test_LQ_lssolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
1605   gsl_test(f, "  LQ_lssolve vander(3)");
1606   s += f;
1607
1608   f = test_LQ_lssolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
1609   gsl_test(f, "  LQ_lssolve vander(4)");
1610   s += f;
1611
1612   f = test_LQ_lssolve_dim(vander12, vander12_solution, 0.05);
1613   gsl_test(f, "  LQ_lssolve vander(12)");
1614   s += f;
1615
1616   return s;
1617 }
1618
1619
1620
1621
1622
1623
1624
1625
1626 int
1627 test_LQ_decomp_dim(const gsl_matrix * m, double eps)
1628 {
1629   int s = 0;
1630   unsigned long i,j, M = m->size1, N = m->size2;
1631
1632   gsl_matrix * lq = gsl_matrix_alloc(M,N);
1633   gsl_matrix * a  = gsl_matrix_alloc(M,N);
1634   gsl_matrix * q  = gsl_matrix_alloc(N,N);
1635   gsl_matrix * l  = gsl_matrix_alloc(M,N);
1636   gsl_vector * d = gsl_vector_alloc(GSL_MIN(M,N));
1637
1638   gsl_matrix_memcpy(lq,m);
1639
1640   s += gsl_linalg_LQ_decomp(lq, d);
1641   s += gsl_linalg_LQ_unpack(lq, d, q, l);
1642
1643    /* compute a = q r */
1644   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, l, q, 0.0, a);
1645
1646   for(i=0; i<M; i++) {
1647     for(j=0; j<N; j++) {
1648       double aij = gsl_matrix_get(a, i, j);
1649       double mij = gsl_matrix_get(m, i, j);
1650       int foo = check(aij, mij, eps);
1651       if(foo) {
1652         printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
1653       }
1654       s += foo;
1655     }
1656   }
1657
1658   gsl_vector_free(d);
1659   gsl_matrix_free(lq);
1660   gsl_matrix_free(a);
1661   gsl_matrix_free(q);
1662   gsl_matrix_free(l);
1663
1664   return s;
1665 }
1666
1667 int test_LQ_decomp(void)
1668 {
1669   int f;
1670   int s = 0;
1671
1672   f = test_LQ_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
1673   gsl_test(f, "  LQ_decomp m(3,5)");
1674   s += f;
1675
1676   f = test_LQ_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
1677   gsl_test(f, "  LQ_decomp m(5,3)");
1678   s += f;
1679
1680   f = test_LQ_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
1681   gsl_test(f, "  LQ_decomp hilbert(2)");
1682   s += f;
1683
1684   f = test_LQ_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
1685   gsl_test(f, "  LQ_decomp hilbert(3)");
1686   s += f;
1687
1688   f = test_LQ_decomp_dim(hilb4, 4 * 1024.0 * GSL_DBL_EPSILON);
1689   gsl_test(f, "  LQ_decomp hilbert(4)");
1690   s += f;
1691
1692   f = test_LQ_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
1693   gsl_test(f, "  LQ_decomp hilbert(12)");
1694   s += f;
1695
1696   f = test_LQ_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
1697   gsl_test(f, "  LQ_decomp vander(2)");
1698   s += f;
1699
1700   f = test_LQ_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
1701   gsl_test(f, "  LQ_decomp vander(3)");
1702   s += f;
1703
1704   f = test_LQ_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
1705   gsl_test(f, "  LQ_decomp vander(4)");
1706   s += f;
1707
1708   f = test_LQ_decomp_dim(vander12, 0.0005);  /* FIXME: bad accuracy */
1709   gsl_test(f, "  LQ_decomp vander(12)");
1710   s += f;
1711
1712   return s;
1713 }
1714
1715
1716
1717
1718 int
1719 test_PTLQ_solve_dim(const gsl_matrix * m, const double * actual, double eps)
1720 {
1721   int s = 0;
1722   int signum;
1723   unsigned long i, dim = m->size1;
1724
1725   gsl_permutation * perm = gsl_permutation_alloc(dim);
1726   gsl_vector * rhs = gsl_vector_alloc(dim);
1727   gsl_matrix * lq  = gsl_matrix_alloc(dim,dim);
1728   gsl_vector * d = gsl_vector_alloc(dim);
1729   gsl_vector * x = gsl_vector_alloc(dim);
1730   gsl_vector * norm = gsl_vector_alloc(dim);
1731
1732   gsl_matrix_transpose_memcpy(lq,m);
1733   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
1734   s += gsl_linalg_PTLQ_decomp(lq, d, perm, &signum, norm);
1735   s += gsl_linalg_PTLQ_solve_T(lq, d, perm, rhs, x);
1736   for(i=0; i<dim; i++) {
1737     int foo = check(gsl_vector_get(x, i), actual[i], eps);
1738     if(foo) {
1739       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
1740     }
1741     s += foo;
1742   }
1743
1744   gsl_vector_free(norm);
1745   gsl_vector_free(x);
1746   gsl_vector_free(d);
1747   gsl_matrix_free(lq);
1748   gsl_vector_free(rhs);
1749   gsl_permutation_free(perm);
1750
1751   return s;
1752 }
1753
1754 int test_PTLQ_solve(void)
1755 {
1756   int f;
1757   int s = 0;
1758
1759   f = test_PTLQ_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
1760   gsl_test(f, "  PTLQ_solve hilbert(2)");
1761   s += f;
1762
1763   f = test_PTLQ_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
1764   gsl_test(f, "  PTLQ_solve hilbert(3)");
1765   s += f;
1766
1767   f = test_PTLQ_solve_dim(hilb4, hilb4_solution, 2 * 2048.0 * GSL_DBL_EPSILON);
1768   gsl_test(f, "  PTLQ_solve hilbert(4)");
1769   s += f;
1770
1771   f = test_PTLQ_solve_dim(hilb12, hilb12_solution, 0.5);
1772   gsl_test(f, "  PTLQ_solve hilbert(12)");
1773   s += f;
1774
1775   f = test_PTLQ_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
1776   gsl_test(f, "  PTLQ_solve vander(2)");
1777   s += f;
1778
1779   f = test_PTLQ_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
1780   gsl_test(f, "  PTLQ_solve vander(3)");
1781   s += f;
1782
1783   f = test_PTLQ_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
1784   gsl_test(f, "  PTLQ_solve vander(4)");
1785   s += f;
1786
1787   f = test_PTLQ_solve_dim(vander12, vander12_solution, 0.05);
1788   gsl_test(f, "  PTLQ_solve vander(12)");
1789   s += f;
1790
1791   return s;
1792 }
1793
1794
1795 int
1796 test_PTLQ_LQsolve_dim(const gsl_matrix * m, const double * actual, double eps)
1797 {
1798   int s = 0;
1799   int signum;
1800   unsigned long i, dim = m->size1;
1801
1802   gsl_permutation * perm = gsl_permutation_alloc(dim);
1803   gsl_vector * rhs = gsl_vector_alloc(dim);
1804   gsl_matrix * lq  = gsl_matrix_alloc(dim,dim);
1805   gsl_matrix * q  = gsl_matrix_alloc(dim,dim);
1806   gsl_matrix * l  = gsl_matrix_alloc(dim,dim);
1807   gsl_vector * d = gsl_vector_alloc(dim);
1808   gsl_vector * x = gsl_vector_alloc(dim);
1809   gsl_vector * norm = gsl_vector_alloc(dim);
1810
1811   gsl_matrix_transpose_memcpy(lq,m);
1812   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
1813   s += gsl_linalg_PTLQ_decomp2(lq, q, l, d, perm, &signum, norm);
1814   s += gsl_linalg_PTLQ_LQsolve_T(q, l, perm, rhs, x);
1815   for(i=0; i<dim; i++) {
1816     int foo = check(gsl_vector_get(x, i), actual[i], eps);
1817     if(foo) {
1818       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
1819     }
1820     s += foo;
1821   }
1822
1823   gsl_vector_free(norm);
1824   gsl_vector_free(x);
1825   gsl_vector_free(d);
1826   gsl_matrix_free(lq);
1827   gsl_vector_free(rhs);
1828   gsl_permutation_free(perm);
1829
1830   return s;
1831 }
1832
1833 int test_PTLQ_LQsolve(void)
1834 {
1835   int f;
1836   int s = 0;
1837
1838   f = test_PTLQ_LQsolve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
1839   gsl_test(f, "  PTLQ_LQsolve hilbert(2)");
1840   s += f;
1841
1842   f = test_PTLQ_LQsolve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
1843   gsl_test(f, "  PTLQ_LQsolve hilbert(3)");
1844   s += f;
1845
1846   f = test_PTLQ_LQsolve_dim(hilb4, hilb4_solution, 2 * 2048.0 * GSL_DBL_EPSILON);
1847   gsl_test(f, "  PTLQ_LQsolve hilbert(4)");
1848   s += f;
1849
1850   f = test_PTLQ_LQsolve_dim(hilb12, hilb12_solution, 0.5);
1851   gsl_test(f, "  PTLQ_LQsolve hilbert(12)");
1852   s += f;
1853
1854   f = test_PTLQ_LQsolve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
1855   gsl_test(f, "  PTLQ_LQsolve vander(2)");
1856   s += f;
1857
1858   f = test_PTLQ_LQsolve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
1859   gsl_test(f, "  PTLQ_LQsolve vander(3)");
1860   s += f;
1861
1862   f = test_PTLQ_LQsolve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
1863   gsl_test(f, "  PTLQ_LQsolve vander(4)");
1864   s += f;
1865
1866   f = test_PTLQ_LQsolve_dim(vander12, vander12_solution, 0.05);
1867   gsl_test(f, "  PTLQ_LQsolve vander(12)");
1868   s += f;
1869
1870   return s;
1871 }
1872
1873
1874 int
1875 test_PTLQ_decomp_dim(const gsl_matrix * m, double eps)
1876 {
1877   int s = 0, signum;
1878   unsigned long i,j, M = m->size1, N = m->size2;
1879
1880   gsl_matrix * lq = gsl_matrix_alloc(N,M);
1881   gsl_matrix * a  = gsl_matrix_alloc(N,M);
1882   gsl_matrix * q  = gsl_matrix_alloc(M,M);
1883   gsl_matrix * l  = gsl_matrix_alloc(N,M);
1884   gsl_vector * d = gsl_vector_alloc(GSL_MIN(M,N));
1885   gsl_vector * norm = gsl_vector_alloc(N);
1886
1887   gsl_permutation * perm = gsl_permutation_alloc(N);
1888
1889   gsl_matrix_transpose_memcpy(lq,m);
1890
1891   s += gsl_linalg_PTLQ_decomp(lq, d, perm, &signum, norm);
1892   s += gsl_linalg_LQ_unpack(lq, d, q, l);
1893
1894    /* compute a = l q */
1895   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, l, q, 0.0, a);
1896
1897
1898    /* Compute P LQ  by permuting the rows of LQ */
1899
1900   for (i = 0; i < M; i++) {
1901     gsl_vector_view col = gsl_matrix_column (a, i);
1902     gsl_permute_vector_inverse (perm, &col.vector);
1903   }
1904
1905   for(i=0; i<M; i++) {
1906     for(j=0; j<N; j++) {
1907       double aij = gsl_matrix_get(a, j, i);
1908       double mij = gsl_matrix_get(m, i, j);
1909       int foo = check(aij, mij, eps);
1910       if(foo) {
1911         printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
1912       }
1913       s += foo;
1914     }
1915   }
1916
1917   gsl_permutation_free (perm);
1918   gsl_vector_free(norm);
1919   gsl_vector_free(d);
1920   gsl_matrix_free(lq);
1921   gsl_matrix_free(a);
1922   gsl_matrix_free(q);
1923   gsl_matrix_free(l);
1924
1925   return s;
1926 }
1927
1928 int test_PTLQ_decomp(void)
1929 {
1930   int f;
1931   int s = 0;
1932
1933   f = test_PTLQ_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
1934   gsl_test(f, "  PTLQ_decomp m(3,5)");
1935   s += f;
1936
1937   f = test_PTLQ_decomp_dim(m53, 2 * 8.0 * GSL_DBL_EPSILON);
1938   gsl_test(f, "  PTLQ_decomp m(5,3)");
1939   s += f;
1940
1941   f = test_PTLQ_decomp_dim(s35, 2 * 8.0 * GSL_DBL_EPSILON);
1942   gsl_test(f, "  PTLQ_decomp s(3,5)");
1943   s += f;
1944
1945   f = test_PTLQ_decomp_dim(s53, 2 * 8.0 * GSL_DBL_EPSILON);
1946   gsl_test(f, "  PTLQ_decomp s(5,3)");
1947   s += f;
1948
1949   f = test_PTLQ_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
1950   gsl_test(f, "  PTLQ_decomp hilbert(2)");
1951   s += f;
1952
1953   f = test_PTLQ_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
1954   gsl_test(f, "  PTLQ_decomp hilbert(3)");
1955   s += f;
1956
1957   f = test_PTLQ_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
1958   gsl_test(f, "  PTLQ_decomp hilbert(4)");
1959   s += f;
1960
1961   f = test_PTLQ_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
1962   gsl_test(f, "  PTLQ_decomp hilbert(12)");
1963   s += f;
1964
1965   f = test_PTLQ_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
1966   gsl_test(f, "  PTLQ_decomp vander(2)");
1967   s += f;
1968
1969   f = test_PTLQ_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
1970   gsl_test(f, "  PTLQ_decomp vander(3)");
1971   s += f;
1972
1973   f = test_PTLQ_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
1974   gsl_test(f, "  PTLQ_decomp vander(4)");
1975   s += f;
1976
1977   f = test_PTLQ_decomp_dim(vander12, 0.0005);  /* FIXME: bad accuracy */
1978   gsl_test(f, "  PTLQ_decomp vander(12)");
1979   s += f;
1980
1981   return s;
1982 }
1983
1984
1985 int
1986 test_LQ_update_dim(const gsl_matrix * m, double eps)
1987 {
1988   int s = 0;
1989   unsigned long i,j, M = m->size1, N = m->size2;
1990
1991   gsl_matrix * lq1  = gsl_matrix_alloc(N,M);
1992   gsl_matrix * lq2  = gsl_matrix_alloc(N,M);
1993   gsl_matrix * q1  = gsl_matrix_alloc(M,M);
1994   gsl_matrix * l1  = gsl_matrix_alloc(N,M);
1995   gsl_matrix * q2  = gsl_matrix_alloc(M,M);
1996   gsl_matrix * l2  = gsl_matrix_alloc(N,M);
1997   gsl_vector * d2 = gsl_vector_alloc(GSL_MIN(M,N));
1998   gsl_vector * u = gsl_vector_alloc(M);
1999   gsl_vector * v = gsl_vector_alloc(N);
2000   gsl_vector * w = gsl_vector_alloc(M);
2001
2002   gsl_matrix_transpose_memcpy(lq1,m);
2003   gsl_matrix_transpose_memcpy(lq2,m);
2004   for(i=0; i<M; i++) gsl_vector_set(u, i, sin(i+1.0));
2005   for(i=0; i<N; i++) gsl_vector_set(v, i, cos(i+2.0) + sin(i*i+3.0));
2006
2007   /* lq1 is updated */
2008
2009   gsl_blas_dger(1.0, v, u, lq1);
2010
2011   /* lq2 is first decomposed, updated later */
2012
2013   s += gsl_linalg_LQ_decomp(lq2, d2);
2014   s += gsl_linalg_LQ_unpack(lq2, d2, q2, l2);
2015
2016   /* compute w = Q^T u */
2017
2018   gsl_blas_dgemv(CblasNoTrans, 1.0, q2, u, 0.0, w);
2019
2020   /* now lq2 is updated */
2021
2022   s += gsl_linalg_LQ_update(q2, l2, v, w);
2023
2024   /* multiply q2*l2 */
2025
2026   gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,l2,q2,0.0,lq2);
2027
2028   /*  check lq1==lq2 */
2029
2030   for(i=0; i<N; i++) {
2031     for(j=0; j<M; j++) {
2032       double s1 = gsl_matrix_get(lq1, i, j);
2033       double s2 = gsl_matrix_get(lq2, i, j);
2034
2035       int foo = check(s1, s2, eps);
2036 #if 0
2037       if(foo) {
2038           printf("LQ:(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, s1, s2);
2039       }
2040 #endif
2041       s += foo;
2042     }
2043   }
2044
2045   gsl_vector_free(d2);
2046   gsl_vector_free(u);
2047   gsl_vector_free(v);
2048   gsl_vector_free(w);
2049   gsl_matrix_free(lq1);
2050   gsl_matrix_free(lq2);
2051   gsl_matrix_free(q1);
2052   gsl_matrix_free(l1);
2053   gsl_matrix_free(q2);
2054   gsl_matrix_free(l2);
2055
2056   return s;
2057 }
2058
2059 int test_LQ_update(void)
2060 {
2061   int f;
2062   int s = 0;
2063
2064   f = test_LQ_update_dim(m35, 2 * 512.0 * GSL_DBL_EPSILON);
2065   gsl_test(f, "  LQ_update m(3,5)");
2066   s += f;
2067
2068   f = test_LQ_update_dim(m53, 2 * 512.0 * GSL_DBL_EPSILON);
2069   gsl_test(f, "  LQ_update m(5,3)");
2070   s += f;
2071
2072   f = test_LQ_update_dim(hilb2,  2 * 512.0 * GSL_DBL_EPSILON);
2073   gsl_test(f, "  LQ_update hilbert(2)");
2074   s += f;
2075
2076   f = test_LQ_update_dim(hilb3,  2 * 512.0 * GSL_DBL_EPSILON);
2077   gsl_test(f, "  LQ_update hilbert(3)");
2078   s += f;
2079
2080   f = test_LQ_update_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
2081   gsl_test(f, "  LQ_update hilbert(4)");
2082   s += f;
2083
2084   f = test_LQ_update_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
2085   gsl_test(f, "  LQ_update hilbert(12)");
2086   s += f;
2087
2088   f = test_LQ_update_dim(vander2, 8.0 * GSL_DBL_EPSILON);
2089   gsl_test(f, "  LQ_update vander(2)");
2090   s += f;
2091
2092   f = test_LQ_update_dim(vander3, 64.0 * GSL_DBL_EPSILON);
2093   gsl_test(f, "  LQ_update vander(3)");
2094   s += f;
2095
2096   f = test_LQ_update_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
2097   gsl_test(f, "  LQ_update vander(4)");
2098   s += f;
2099
2100   f = test_LQ_update_dim(vander12, 0.0005);  /* FIXME: bad accuracy */
2101   gsl_test(f, "  LQ_update vander(12)");
2102   s += f;
2103
2104   return s;
2105 }
2106
2107 int
2108 test_SV_solve_dim(const gsl_matrix * m, const double * actual, double eps)
2109 {
2110   int s = 0;
2111   unsigned long i, dim = m->size1;
2112
2113   gsl_vector * rhs = gsl_vector_alloc(dim);
2114   gsl_matrix * u  = gsl_matrix_alloc(dim,dim);
2115   gsl_matrix * q  = gsl_matrix_alloc(dim,dim);
2116   gsl_vector * d = gsl_vector_alloc(dim);
2117   gsl_vector * x = gsl_vector_calloc(dim);
2118   gsl_matrix_memcpy(u,m);
2119   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
2120   s += gsl_linalg_SV_decomp(u, q, d, x);
2121   s += gsl_linalg_SV_solve(u, q, d, rhs, x);
2122   for(i=0; i<dim; i++) {
2123     int foo = check(gsl_vector_get(x, i), actual[i], eps);
2124     if(foo) {
2125       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
2126     }
2127     s += foo;
2128   }
2129   gsl_vector_free(x);
2130   gsl_vector_free(d);
2131   gsl_matrix_free(u);
2132   gsl_matrix_free(q);
2133   gsl_vector_free(rhs);
2134
2135   return s;
2136 }
2137
2138 int test_SV_solve(void)
2139 {
2140   int f;
2141   int s = 0;
2142
2143   f = test_SV_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
2144   gsl_test(f, "  SV_solve hilbert(2)");
2145   s += f;
2146
2147   f = test_SV_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
2148   gsl_test(f, "  SV_solve hilbert(3)");
2149   s += f;
2150
2151   f = test_SV_solve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
2152   gsl_test(f, "  SV_solve hilbert(4)");
2153   s += f;
2154
2155   f = test_SV_solve_dim(hilb12, hilb12_solution, 0.5);
2156   gsl_test(f, "  SV_solve hilbert(12)");
2157   s += f;
2158
2159   f = test_SV_solve_dim(vander2, vander2_solution, 64.0 * GSL_DBL_EPSILON);
2160   gsl_test(f, "  SV_solve vander(2)");
2161   s += f;
2162
2163   f = test_SV_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
2164   gsl_test(f, "  SV_solve vander(3)");
2165   s += f;
2166
2167   f = test_SV_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
2168   gsl_test(f, "  SV_solve vander(4)");
2169   s += f;
2170
2171   f = test_SV_solve_dim(vander12, vander12_solution, 0.05);
2172   gsl_test(f, "  SV_solve vander(12)");
2173   s += f;
2174
2175   return s;
2176 }
2177
2178 int
2179 test_SV_decomp_dim(const gsl_matrix * m, double eps)
2180 {
2181   int s = 0;
2182   double di1;
2183   unsigned long i,j, M = m->size1, N = m->size2;
2184
2185   gsl_matrix * v  = gsl_matrix_alloc(M,N);
2186   gsl_matrix * a  = gsl_matrix_alloc(M,N);
2187   gsl_matrix * q  = gsl_matrix_alloc(N,N);
2188   gsl_matrix * dqt  = gsl_matrix_alloc(N,N);
2189   gsl_vector * d  = gsl_vector_alloc(N);
2190   gsl_vector * w  = gsl_vector_alloc(N);
2191
2192   gsl_matrix_memcpy(v,m);
2193
2194   s += gsl_linalg_SV_decomp(v, q, d, w);
2195
2196   /* Check that singular values are non-negative and in non-decreasing
2197      order */
2198
2199   di1 = 0.0;
2200
2201   for (i = 0; i < N; i++)
2202     {
2203       double di = gsl_vector_get (d, i);
2204
2205       if (gsl_isnan (di))
2206         {
2207           continue;  /* skip NaNs */
2208         }
2209
2210       if (di < 0) {
2211         s++;
2212         printf("singular value %lu = %22.18g < 0\n", i, di);
2213       }
2214
2215       if(i > 0 && di > di1) {
2216         s++;
2217         printf("singular value %lu = %22.18g vs previous %22.18g\n", i, di, di1);
2218       }
2219
2220       di1 = di;
2221     }
2222
2223   /* Scale dqt = D Q^T */
2224
2225   for (i = 0; i < N ; i++)
2226     {
2227       double di = gsl_vector_get (d, i);
2228
2229       for (j = 0; j < N; j++)
2230         {
2231           double qji = gsl_matrix_get(q, j, i);
2232           gsl_matrix_set (dqt, i, j, qji * di);
2233         }
2234     }
2235
2236   /* compute a = v dqt */
2237   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, v, dqt, 0.0, a);
2238
2239   for(i=0; i<M; i++) {
2240     for(j=0; j<N; j++) {
2241       double aij = gsl_matrix_get(a, i, j);
2242       double mij = gsl_matrix_get(m, i, j);
2243       int foo = check(aij, mij, eps);
2244       if(foo) {
2245         printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
2246       }
2247       s += foo;
2248     }
2249   }
2250   gsl_vector_free(w);
2251   gsl_vector_free(d);
2252   gsl_matrix_free(v);
2253   gsl_matrix_free(a);
2254   gsl_matrix_free(q);
2255   gsl_matrix_free(dqt);
2256
2257   return s;
2258 }
2259
2260 int test_SV_decomp(void)
2261 {
2262   int f;
2263   int s = 0;
2264
2265   f = test_SV_decomp_dim(m11, 2 * GSL_DBL_EPSILON);
2266   gsl_test(f, "  SV_decomp m(1,1)");
2267   s += f;
2268
2269   f = test_SV_decomp_dim(m51, 2 * 64.0 * GSL_DBL_EPSILON);
2270   gsl_test(f, "  SV_decomp m(5,1)");
2271   s += f;
2272
2273   /* M<N not implemented yet */
2274 #if 0
2275   f = test_SV_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
2276   gsl_test(f, "  SV_decomp m(3,5)");
2277   s += f;
2278 #endif
2279   f = test_SV_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
2280   gsl_test(f, "  SV_decomp m(5,3)");
2281   s += f;
2282
2283   f = test_SV_decomp_dim(moler10, 2 * 64.0 * GSL_DBL_EPSILON);
2284   gsl_test(f, "  SV_decomp moler(10)");
2285   s += f;
2286
2287   f = test_SV_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
2288   gsl_test(f, "  SV_decomp hilbert(2)");
2289   s += f;
2290
2291   f = test_SV_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
2292   gsl_test(f, "  SV_decomp hilbert(3)");
2293   s += f;
2294
2295   f = test_SV_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
2296   gsl_test(f, "  SV_decomp hilbert(4)");
2297   s += f;
2298
2299   f = test_SV_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
2300   gsl_test(f, "  SV_decomp hilbert(12)");
2301   s += f;
2302
2303   f = test_SV_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
2304   gsl_test(f, "  SV_decomp vander(2)");
2305   s += f;
2306
2307   f = test_SV_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
2308   gsl_test(f, "  SV_decomp vander(3)");
2309   s += f;
2310
2311   f = test_SV_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
2312   gsl_test(f, "  SV_decomp vander(4)");
2313   s += f;
2314
2315   f = test_SV_decomp_dim(vander12, 1e-4);
2316   gsl_test(f, "  SV_decomp vander(12)");
2317   s += f;
2318
2319   f = test_SV_decomp_dim(row3, 10 * GSL_DBL_EPSILON);
2320   gsl_test(f, "  SV_decomp row3");
2321   s += f;
2322
2323   f = test_SV_decomp_dim(row5, 128 * GSL_DBL_EPSILON);
2324   gsl_test(f, "  SV_decomp row5");
2325   s += f;
2326
2327   f = test_SV_decomp_dim(row12, 1024 * GSL_DBL_EPSILON);
2328   gsl_test(f, "  SV_decomp row12");
2329   s += f;
2330
2331   f = test_SV_decomp_dim(inf5, 1024 * GSL_DBL_EPSILON);
2332   gsl_test(f, "  SV_decomp inf5");
2333   s += f;
2334
2335   f = test_SV_decomp_dim(nan5, 1024 * GSL_DBL_EPSILON);
2336   gsl_test(f, "  SV_decomp nan5");
2337   s += f;
2338
2339   f = test_SV_decomp_dim(dblmin3, 1024 * GSL_DBL_EPSILON);
2340   gsl_test(f, "  SV_decomp dblmin3");
2341   s += f;
2342
2343   f = test_SV_decomp_dim(dblmin5, 1024 * GSL_DBL_EPSILON);
2344   gsl_test(f, "  SV_decomp dblmin5");
2345   s += f;
2346
2347
2348   {
2349     double i1, i2, i3, i4;
2350     double lower = -2, upper = 2;
2351
2352     for (i1 = lower; i1 <= upper; i1++)
2353       {
2354         for (i2 = lower; i2 <= upper; i2++)
2355           {
2356             for (i3 = lower; i3 <= upper; i3++)
2357               {
2358                 for (i4 = lower; i4 <= upper; i4++)
2359                   {
2360                     gsl_matrix_set (A22, 0,0, i1);
2361                     gsl_matrix_set (A22, 0,1, i2);
2362                     gsl_matrix_set (A22, 1,0, i3);
2363                     gsl_matrix_set (A22, 1,1, i4);
2364
2365                     f = test_SV_decomp_dim(A22, 16 * GSL_DBL_EPSILON);
2366                     gsl_test(f, "  SV_decomp (2x2) A=[%g, %g; %g, %g]", i1,i2,i3,i4);
2367                     s += f;
2368                   }
2369               }
2370           }
2371       }
2372   }
2373
2374   {
2375     int i;
2376     double carry = 0, lower = 0, upper = 1;
2377     double *a = A33->data;
2378
2379     for (i=0; i<9; i++) {
2380       a[i] = lower;
2381     }
2382
2383     while (carry == 0.0) {
2384       f = test_SV_decomp_dim(A33, 64 * GSL_DBL_EPSILON);
2385       gsl_test(f, "  SV_decomp (3x3) A=[ %g, %g, %g; %g, %g, %g; %g, %g, %g]",
2386                a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8]);
2387
2388       /* increment */
2389       carry=1.0;
2390       for (i=9; carry > 0.0 && i>0 && i--;)
2391         {
2392           double v=a[i]+carry;
2393           carry = (v>upper) ? 1.0 : 0.0;
2394           a[i] = (v>upper) ? lower : v;
2395         }
2396     }
2397   }
2398
2399 #ifdef TEST_SVD_4X4
2400   {
2401     int i;
2402     double carry = 0, lower = 0, upper = 1;
2403     double *a = A44->data;
2404
2405     for (i=0; i<16; i++) {
2406       a[i] = lower;
2407     }
2408
2409     while (carry == 0.0) {
2410       f = test_SV_decomp_dim(A44, 64 * GSL_DBL_EPSILON);
2411       gsl_test(f, "  SV_decomp (4x4) A=[ %g, %g, %g, %g; %g, %g, %g, %g; %g, %g, %g, %g; %g, %g, %g, %g]",
2412                a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9],
2413                a[10], a[11], a[12], a[13], a[14], a[15]);
2414
2415       /* increment */
2416       carry=1.0;
2417       for (i=16; carry > 0.0 && i>0 && i--;)
2418         {
2419           double v=a[i]+carry;
2420           carry = (v>upper) ? 1.0 : 0.0;
2421           a[i] = (v>upper) ? lower : v;
2422         }
2423     }
2424   }
2425 #endif
2426
2427   return s;
2428 }
2429
2430
2431 int
2432 test_SV_decomp_mod_dim(const gsl_matrix * m, double eps)
2433 {
2434   int s = 0;
2435   double di1;
2436   unsigned long i,j, M = m->size1, N = m->size2;
2437
2438   gsl_matrix * v  = gsl_matrix_alloc(M,N);
2439   gsl_matrix * a  = gsl_matrix_alloc(M,N);
2440   gsl_matrix * q  = gsl_matrix_alloc(N,N);
2441   gsl_matrix * x  = gsl_matrix_alloc(N,N);
2442   gsl_matrix * dqt  = gsl_matrix_alloc(N,N);
2443   gsl_vector * d  = gsl_vector_alloc(N);
2444   gsl_vector * w  = gsl_vector_alloc(N);
2445
2446   gsl_matrix_memcpy(v,m);
2447
2448   s += gsl_linalg_SV_decomp_mod(v, x, q, d, w);
2449
2450   /* Check that singular values are non-negative and in non-decreasing
2451      order */
2452
2453   di1 = 0.0;
2454
2455   for (i = 0; i < N; i++)
2456     {
2457       double di = gsl_vector_get (d, i);
2458
2459       if (gsl_isnan (di))
2460         {
2461           continue;  /* skip NaNs */
2462         }
2463
2464       if (di < 0) {
2465         s++;
2466         printf("singular value %lu = %22.18g < 0\n", i, di);
2467       }
2468
2469       if(i > 0 && di > di1) {
2470         s++;
2471         printf("singular value %lu = %22.18g vs previous %22.18g\n", i, di, di1);
2472       }
2473
2474       di1 = di;
2475     }
2476
2477   /* Scale dqt = D Q^T */
2478
2479   for (i = 0; i < N ; i++)
2480     {
2481       double di = gsl_vector_get (d, i);
2482
2483       for (j = 0; j < N; j++)
2484         {
2485           double qji = gsl_matrix_get(q, j, i);
2486           gsl_matrix_set (dqt, i, j, qji * di);
2487         }
2488     }
2489
2490   /* compute a = v dqt */
2491   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, v, dqt, 0.0, a);
2492
2493   for(i=0; i<M; i++) {
2494     for(j=0; j<N; j++) {
2495       double aij = gsl_matrix_get(a, i, j);
2496       double mij = gsl_matrix_get(m, i, j);
2497       int foo = check(aij, mij, eps);
2498       if(foo) {
2499         printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
2500       }
2501       s += foo;
2502     }
2503   }
2504   gsl_vector_free(w);
2505   gsl_vector_free(d);
2506   gsl_matrix_free(v);
2507   gsl_matrix_free(a);
2508   gsl_matrix_free(q);
2509   gsl_matrix_free(dqt);
2510   gsl_matrix_free (x);
2511
2512   return s;
2513 }
2514
2515 int test_SV_decomp_mod(void)
2516 {
2517   int f;
2518   int s = 0;
2519
2520   f = test_SV_decomp_mod_dim(m11, 2 * GSL_DBL_EPSILON);
2521   gsl_test(f, "  SV_decomp_mod m(1,1)");
2522   s += f;
2523
2524   f = test_SV_decomp_mod_dim(m51, 2 * 64.0 * GSL_DBL_EPSILON);
2525   gsl_test(f, "  SV_decomp_mod m(5,1)");
2526   s += f;
2527
2528   /* M<N not implemented yet */
2529 #if 0
2530   f = test_SV_decomp_mod_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
2531   gsl_test(f, "  SV_decomp_mod m(3,5)");
2532   s += f;
2533 #endif
2534   f = test_SV_decomp_mod_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
2535   gsl_test(f, "  SV_decomp_mod m(5,3)");
2536   s += f;
2537
2538   f = test_SV_decomp_mod_dim(moler10, 2 * 64.0 * GSL_DBL_EPSILON);
2539   gsl_test(f, "  SV_decomp_mod moler(10)");
2540   s += f;
2541
2542   f = test_SV_decomp_mod_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
2543   gsl_test(f, "  SV_decomp_mod hilbert(2)");
2544   s += f;
2545
2546   f = test_SV_decomp_mod_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
2547   gsl_test(f, "  SV_decomp_mod hilbert(3)");
2548   s += f;
2549
2550   f = test_SV_decomp_mod_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
2551   gsl_test(f, "  SV_decomp_mod hilbert(4)");
2552   s += f;
2553
2554   f = test_SV_decomp_mod_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
2555   gsl_test(f, "  SV_decomp_mod hilbert(12)");
2556   s += f;
2557
2558   f = test_SV_decomp_mod_dim(vander2, 8.0 * GSL_DBL_EPSILON);
2559   gsl_test(f, "  SV_decomp_mod vander(2)");
2560   s += f;
2561
2562   f = test_SV_decomp_mod_dim(vander3, 64.0 * GSL_DBL_EPSILON);
2563   gsl_test(f, "  SV_decomp_mod vander(3)");
2564   s += f;
2565
2566   f = test_SV_decomp_mod_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
2567   gsl_test(f, "  SV_decomp_mod vander(4)");
2568   s += f;
2569
2570   f = test_SV_decomp_mod_dim(vander12, 1e-4);
2571   gsl_test(f, "  SV_decomp_mod vander(12)");
2572   s += f;
2573
2574   f = test_SV_decomp_mod_dim(row3, 10 * GSL_DBL_EPSILON);
2575   gsl_test(f, "  SV_decomp_mod row3");
2576   s += f;
2577
2578   f = test_SV_decomp_mod_dim(row5, 128 * GSL_DBL_EPSILON);
2579   gsl_test(f, "  SV_decomp_mod row5");
2580   s += f;
2581
2582   f = test_SV_decomp_mod_dim(row12, 1024 * GSL_DBL_EPSILON);
2583   gsl_test(f, "  SV_decomp_mod row12");
2584   s += f;
2585
2586   f = test_SV_decomp_mod_dim(inf5, 1024 * GSL_DBL_EPSILON);
2587   gsl_test(f, "  SV_decomp_mod inf5");
2588   s += f;
2589
2590   f = test_SV_decomp_mod_dim(nan5, 1024 * GSL_DBL_EPSILON);
2591   gsl_test(f, "  SV_decomp_mod nan5");
2592   s += f;
2593
2594
2595   {
2596     double i1, i2, i3, i4;
2597     double lower = -2, upper = 2;
2598
2599     for (i1 = lower; i1 <= upper; i1++)
2600       {
2601         for (i2 = lower; i2 <= upper; i2++)
2602           {
2603             for (i3 = lower; i3 <= upper; i3++)
2604               {
2605                 for (i4 = lower; i4 <= upper; i4++)
2606                   {
2607                     gsl_matrix_set (A22, 0,0, i1);
2608                     gsl_matrix_set (A22, 0,1, i2);
2609                     gsl_matrix_set (A22, 1,0, i3);
2610                     gsl_matrix_set (A22, 1,1, i4);
2611
2612                     f = test_SV_decomp_mod_dim(A22, 16 * GSL_DBL_EPSILON);
2613                     gsl_test(f, "  SV_decomp_mod (2x2) A=[%g, %g; %g, %g]", i1,i2,i3,i4);
2614                     s += f;
2615                   }
2616               }
2617           }
2618       }
2619   }
2620
2621   {
2622     int i;
2623     double carry = 0, lower = 0, upper = 1;
2624     double *a = A33->data;
2625
2626     for (i=0; i<9; i++) {
2627       a[i] = lower;
2628     }
2629
2630     while (carry == 0.0) {
2631       f = test_SV_decomp_mod_dim(A33, 64 * GSL_DBL_EPSILON);
2632       gsl_test(f, "  SV_decomp_mod (3x3) A=[ %g, %g, %g; %g, %g, %g; %g, %g, %g]",
2633                a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8]);
2634
2635       /* increment */
2636       carry=1.0;
2637       for (i=9; carry > 0.0 && i>0 && i--;)
2638         {
2639           double v=a[i]+carry;
2640           carry = (v>upper) ? 1.0 : 0.0;
2641           a[i] = (v>upper) ? lower : v;
2642         }
2643     }
2644   }
2645
2646 #ifdef TEST_SVD_4X4
2647   {
2648     int i;
2649     double carry = 0, lower = 0, upper = 1;
2650     double *a = A44->data;
2651
2652     for (i=0; i<16; i++) {
2653       a[i] = lower;
2654     }
2655
2656     while (carry == 0.0) {
2657       f = test_SV_decomp_mod_dim(A44, 64 * GSL_DBL_EPSILON);
2658       gsl_test(f, "  SV_decomp_mod (4x4) A=[ %g, %g, %g, %g; %g, %g, %g, %g; %g, %g, %g, %g; %g, %g, %g, %g]",
2659                a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9],
2660                a[10], a[11], a[12], a[13], a[14], a[15]);
2661
2662       /* increment */
2663       carry=1.0;
2664       for (i=16; carry>0.0 && i>0 && i--;)
2665         {
2666           double v=a[i]+carry;
2667           carry = (v>upper) ? 1.0 : 0.0;
2668           a[i] = (v>upper) ? lower : v;
2669         }
2670     }
2671   }
2672 #endif
2673
2674   return s;
2675 }
2676
2677
2678 int
2679 test_SV_decomp_jacobi_dim(const gsl_matrix * m, double eps)
2680 {
2681   int s = 0;
2682   double di1;
2683   unsigned long i,j, M = m->size1, N = m->size2;
2684
2685   gsl_matrix * v  = gsl_matrix_alloc(M,N);
2686   gsl_matrix * a  = gsl_matrix_alloc(M,N);
2687   gsl_matrix * q  = gsl_matrix_alloc(N,N);
2688   gsl_matrix * dqt  = gsl_matrix_alloc(N,N);
2689   gsl_vector * d  = gsl_vector_alloc(N);
2690
2691   gsl_matrix_memcpy(v,m);
2692
2693   s += gsl_linalg_SV_decomp_jacobi(v, q, d);
2694   if (s)
2695     printf("call returned status = %d\n", s);
2696
2697   /* Check that singular values are non-negative and in non-decreasing
2698      order */
2699
2700   di1 = 0.0;
2701
2702   for (i = 0; i < N; i++)
2703     {
2704       double di = gsl_vector_get (d, i);
2705
2706       if (gsl_isnan (di))
2707         {
2708           continue;  /* skip NaNs */
2709         }
2710
2711       if (di < 0) {
2712         s++;
2713         printf("singular value %lu = %22.18g < 0\n", i, di);
2714       }
2715
2716       if(i > 0 && di > di1) {
2717         s++;
2718         printf("singular value %lu = %22.18g vs previous %22.18g\n", i, di, di1);
2719       }
2720
2721       di1 = di;
2722     }
2723
2724   /* Scale dqt = D Q^T */
2725
2726   for (i = 0; i < N ; i++)
2727     {
2728       double di = gsl_vector_get (d, i);
2729
2730       for (j = 0; j < N; j++)
2731         {
2732           double qji = gsl_matrix_get(q, j, i);
2733           gsl_matrix_set (dqt, i, j, qji * di);
2734         }
2735     }
2736
2737   /* compute a = v dqt */
2738   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, v, dqt, 0.0, a);
2739
2740   for(i=0; i<M; i++) {
2741     for(j=0; j<N; j++) {
2742       double aij = gsl_matrix_get(a, i, j);
2743       double mij = gsl_matrix_get(m, i, j);
2744       int foo = check(aij, mij, eps);
2745       if(foo) {
2746         printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
2747       }
2748       s += foo;
2749     }
2750   }
2751   gsl_vector_free(d);
2752   gsl_matrix_free(v);
2753   gsl_matrix_free(a);
2754   gsl_matrix_free(q);
2755   gsl_matrix_free(dqt);
2756
2757   return s;
2758 }
2759
2760 int test_SV_decomp_jacobi(void)
2761 {
2762   int f;
2763   int s = 0;
2764
2765   f = test_SV_decomp_jacobi_dim(m11, 2 * GSL_DBL_EPSILON);
2766   gsl_test(f, "  SV_decomp_jacobi m(1,1)");
2767   s += f;
2768
2769   f = test_SV_decomp_jacobi_dim(m51, 2 * 64.0 * GSL_DBL_EPSILON);
2770   gsl_test(f, "  SV_decomp_jacobi m(5,1)");
2771   s += f;
2772
2773   /* M<N not implemented yet */
2774 #if 0
2775   f = test_SV_decomp_jacobi_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
2776   gsl_test(f, "  SV_decomp_jacobi m(3,5)");
2777   s += f;
2778 #endif
2779   f = test_SV_decomp_jacobi_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
2780   gsl_test(f, "  SV_decomp_jacobi m(5,3)");
2781   s += f;
2782
2783   f = test_SV_decomp_jacobi_dim(moler10, 2 * 64.0 * GSL_DBL_EPSILON);
2784   gsl_test(f, "  SV_decomp_jacobi moler(10)");
2785   s += f;
2786
2787   f = test_SV_decomp_jacobi_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
2788   gsl_test(f, "  SV_decomp_jacobi hilbert(2)");
2789   s += f;
2790
2791   f = test_SV_decomp_jacobi_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
2792   gsl_test(f, "  SV_decomp_jacobi hilbert(3)");
2793   s += f;
2794
2795   f = test_SV_decomp_jacobi_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
2796   gsl_test(f, "  SV_decomp_jacobi hilbert(4)");
2797   s += f;
2798
2799   f = test_SV_decomp_jacobi_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
2800   gsl_test(f, "  SV_decomp_jacobi hilbert(12)");
2801   s += f;
2802
2803   f = test_SV_decomp_jacobi_dim(vander2, 8.0 * GSL_DBL_EPSILON);
2804   gsl_test(f, "  SV_decomp_jacobi vander(2)");
2805   s += f;
2806
2807   f = test_SV_decomp_jacobi_dim(vander3, 64.0 * GSL_DBL_EPSILON);
2808   gsl_test(f, "  SV_decomp_jacobi vander(3)");
2809   s += f;
2810
2811   f = test_SV_decomp_jacobi_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
2812   gsl_test(f, "  SV_decomp_jacobi vander(4)");
2813   s += f;
2814
2815   f = test_SV_decomp_jacobi_dim(vander12, 1e-4);
2816   gsl_test(f, "  SV_decomp_jacobi vander(12)");
2817   s += f;
2818
2819   f = test_SV_decomp_jacobi_dim(row3, 10 * GSL_DBL_EPSILON);
2820   gsl_test(f, "  SV_decomp_jacobi row3");
2821   s += f;
2822
2823   f = test_SV_decomp_jacobi_dim(row5, 128 * GSL_DBL_EPSILON);
2824   gsl_test(f, "  SV_decomp_jacobi row5");
2825   s += f;
2826
2827   f = test_SV_decomp_jacobi_dim(row12, 1024 * GSL_DBL_EPSILON);
2828   gsl_test(f, "  SV_decomp_jacobi row12");
2829   s += f;
2830
2831
2832 #ifdef TEST_JACOBI_INF
2833   f = test_SV_decomp_jacobi_dim(inf5, 1024 * GSL_DBL_EPSILON);
2834   gsl_test(f, "  SV_decomp_jacobi inf5");
2835   s += f;
2836
2837   f = test_SV_decomp_jacobi_dim(nan5, 1024 * GSL_DBL_EPSILON);
2838   gsl_test(f, "  SV_decomp_jacobi nan5");
2839   s += f;
2840 #endif
2841
2842   {
2843     double i1, i2, i3, i4;
2844     double lower = -2, upper = 2;
2845
2846     for (i1 = lower; i1 <= upper; i1++)
2847       {
2848         for (i2 = lower; i2 <= upper; i2++)
2849           {
2850             for (i3 = lower; i3 <= upper; i3++)
2851               {
2852                 for (i4 = lower; i4 <= upper; i4++)
2853                   {
2854                     gsl_matrix_set (A22, 0,0, i1);
2855                     gsl_matrix_set (A22, 0,1, i2);
2856                     gsl_matrix_set (A22, 1,0, i3);
2857                     gsl_matrix_set (A22, 1,1, i4);
2858
2859                     f = test_SV_decomp_jacobi_dim(A22, 16 * GSL_DBL_EPSILON);
2860                     gsl_test(f, "  SV_decomp_jacobi (2x2) A=[%g, %g; %g, %g]", i1,i2,i3,i4);
2861                     s += f;
2862                   }
2863               }
2864           }
2865       }
2866   }
2867
2868   {
2869     int i;
2870     double carry = 0, lower = 0, upper = 1;
2871     double *a = A33->data;
2872
2873     for (i=0; i<9; i++) {
2874       a[i] = lower;
2875     }
2876
2877     while (carry == 0.0) {
2878       f = test_SV_decomp_jacobi_dim(A33, 64 * GSL_DBL_EPSILON);
2879       gsl_test(f, "  SV_decomp_jacobi (3x3) A=[ %g, %g, %g; %g, %g, %g; %g, %g, %g]",
2880                a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8]);
2881
2882       /* increment */
2883       carry=1.0;
2884       for (i=9; carry > 0.0 && i>0 && i--;)
2885         {
2886           double v=a[i]+carry;
2887           carry = (v>upper) ? 1.0 : 0.0;
2888           a[i] = (v>upper) ? lower : v;
2889         }
2890     }
2891   }
2892
2893 #ifdef TEST_SVD_4X4
2894   {
2895     int i;
2896     unsigned long k = 0;
2897     double carry = 0, lower = 0, upper = 1;
2898     double *a = A44->data;
2899
2900     for (i=0; i<16; i++) {
2901       a[i] = lower;
2902     }
2903
2904     while (carry == 0.0) {
2905       k++;
2906       f = test_SV_decomp_jacobi_dim(A44, 64 * GSL_DBL_EPSILON);
2907       gsl_test(f, "  SV_decomp_jacobi (4x4) A=[ %g, %g, %g, %g; %g, %g, %g, %g; %g, %g, %g, %g; %g, %g, %g, %g] %lu",
2908                a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9],
2909                a[10], a[11], a[12], a[13], a[14], a[15], k);
2910       /* increment */
2911       carry=1.0;
2912       for (i=16; carry > 0.0 && i>0 && i--;)
2913         {
2914           double v=a[i]+carry;
2915           carry = (v>upper) ? 1.0 : 0.0;
2916           a[i] = (v>upper) ? lower : v;
2917         }
2918     }
2919   }
2920 #endif
2921
2922   {
2923     int i;
2924     unsigned long k = 0;
2925     double carry = 0, lower = 0, upper = 1;
2926     double *a = A55->data;
2927
2928     for (i=0; i<25; i++) {
2929       a[i] = lower;
2930     }
2931
2932     while (carry == 0.0) {
2933       k++;
2934
2935       if (k % 1001 == 0)
2936         {
2937           f = test_SV_decomp_jacobi_dim(A55, 64 * GSL_DBL_EPSILON);
2938           gsl_test(f, "  SV_decomp_jacobi (5x5) case=%lu",k);
2939         }
2940
2941       /* increment */
2942       carry=1.0;
2943       for (i=25; carry >0.0 && i>0 && i--;)
2944         {
2945           double v=a[i]+carry;
2946           carry = (v>upper) ? 1.0 : 0.0;
2947           a[i] = (v>upper) ? lower : v;
2948         }
2949     }
2950   }
2951
2952
2953   return s;
2954 }
2955
2956
2957 int
2958 test_cholesky_solve_dim(const gsl_matrix * m, const double * actual, double eps)
2959 {
2960   int s = 0;
2961   unsigned long i, dim = m->size1;
2962
2963   gsl_vector * rhs = gsl_vector_alloc(dim);
2964   gsl_matrix * u  = gsl_matrix_alloc(dim,dim);
2965   gsl_vector * x = gsl_vector_calloc(dim);
2966   gsl_matrix_memcpy(u,m);
2967   for(i=0; i<dim; i++) gsl_vector_set(rhs, i, i+1.0);
2968   s += gsl_linalg_cholesky_decomp(u);
2969   s += gsl_linalg_cholesky_solve(u, rhs, x);
2970   for(i=0; i<dim; i++) {
2971     int foo = check(gsl_vector_get(x, i), actual[i], eps);
2972     if(foo) {
2973       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
2974     }
2975     s += foo;
2976   }
2977   gsl_vector_free(x);
2978   gsl_matrix_free(u);
2979   gsl_vector_free(rhs);
2980
2981   return s;
2982 }
2983
2984 int test_cholesky_solve(void)
2985 {
2986   int f;
2987   int s = 0;
2988
2989   f = test_cholesky_solve_dim(hilb2, hilb2_solution, 2 * 8.0 * GSL_DBL_EPSILON);
2990   gsl_test(f, "  cholesky_solve hilbert(2)");
2991   s += f;
2992
2993   f = test_cholesky_solve_dim(hilb3, hilb3_solution, 2 * 64.0 * GSL_DBL_EPSILON);
2994   gsl_test(f, "  cholesky_solve hilbert(3)");
2995   s += f;
2996
2997   f = test_cholesky_solve_dim(hilb4, hilb4_solution, 2 * 1024.0 * GSL_DBL_EPSILON);
2998   gsl_test(f, "  cholesky_solve hilbert(4)");
2999   s += f;
3000
3001   f = test_cholesky_solve_dim(hilb12, hilb12_solution, 0.5);
3002   gsl_test(f, "  cholesky_solve hilbert(12)");
3003   s += f;
3004
3005   return s;
3006 }
3007
3008
3009 int
3010 test_cholesky_decomp_dim(const gsl_matrix * m, double eps)
3011 {
3012   int s = 0;
3013   unsigned long i,j, M = m->size1, N = m->size2;
3014
3015   gsl_matrix * v  = gsl_matrix_alloc(M,N);
3016   gsl_matrix * a  = gsl_matrix_alloc(M,N);
3017   gsl_matrix * l  = gsl_matrix_alloc(M,N);
3018   gsl_matrix * lt  = gsl_matrix_alloc(N,N);
3019
3020   gsl_matrix_memcpy(v,m);
3021
3022   s += gsl_linalg_cholesky_decomp(v);
3023
3024   /* Compute L LT */
3025
3026   for (i = 0; i < N ; i++)
3027     {
3028       for (j = 0; j < N; j++)
3029         {
3030           double vij = gsl_matrix_get(v, i, j);
3031           gsl_matrix_set (l, i, j, i>=j ? vij : 0);
3032           gsl_matrix_set (lt, i, j, i<=j ? vij : 0);
3033         }
3034     }
3035
3036   /* compute a = l lt */
3037   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, l, lt, 0.0, a);
3038
3039   for(i=0; i<M; i++) {
3040     for(j=0; j<N; j++) {
3041       double aij = gsl_matrix_get(a, i, j);
3042       double mij = gsl_matrix_get(m, i, j);
3043       int foo = check(aij, mij, eps);
3044       if(foo) {
3045         printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
3046       }
3047       s += foo;
3048     }
3049   }
3050
3051   gsl_matrix_free(v);
3052   gsl_matrix_free(a);
3053   gsl_matrix_free(l);
3054   gsl_matrix_free(lt);
3055
3056   return s;
3057 }
3058
3059 int test_cholesky_decomp(void)
3060 {
3061   int f;
3062   int s = 0;
3063
3064   f = test_cholesky_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
3065   gsl_test(f, "  cholesky_decomp hilbert(2)");
3066   s += f;
3067
3068   f = test_cholesky_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
3069   gsl_test(f, "  cholesky_decomp hilbert(3)");
3070   s += f;
3071
3072   f = test_cholesky_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
3073   gsl_test(f, "  cholesky_decomp hilbert(4)");
3074   s += f;
3075
3076   f = test_cholesky_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
3077   gsl_test(f, "  cholesky_decomp hilbert(12)");
3078   s += f;
3079
3080   return s;
3081 }
3082
3083
3084 int
3085 test_cholesky_decomp_unit_dim(const gsl_matrix * m, double eps)
3086 {
3087   int s = 0;
3088   const unsigned long M = m->size1;
3089   const unsigned long N = m->size2;
3090   unsigned long i,j;
3091
3092   gsl_matrix * v  = gsl_matrix_alloc(M,N);
3093   gsl_matrix * a  = gsl_matrix_alloc(M,N);
3094   gsl_matrix * l  = gsl_matrix_alloc(M,N);
3095   gsl_matrix * lt = gsl_matrix_alloc(N,N);
3096   gsl_matrix * dm = gsl_matrix_alloc(M,N);
3097   gsl_vector * dv = gsl_vector_alloc(M);
3098
3099   gsl_matrix_memcpy(v,m);
3100
3101   s += gsl_linalg_cholesky_decomp_unit(v, dv);
3102
3103   /*
3104   for(i = 0; i < M; i++)
3105   {
3106     for(j = 0; j < N; j++)
3107     {
3108       printf("v[%lu,%lu]: %22.18e\n", i,j, gsl_matrix_get(v, i, j));
3109     }
3110   }
3111
3112
3113   for(i = 0; i < M; i++)
3114   {
3115     printf("d[%lu]: %22.18e\n", i, gsl_vector_get(dv, i));
3116   }
3117   */
3118
3119   /* put L and transpose(L) into separate matrices */
3120
3121   for(i = 0; i < N ; i++)
3122   {
3123     for(j = 0; j < N; j++)
3124     {
3125       const double vij = gsl_matrix_get(v, i, j);
3126       gsl_matrix_set (l,  i, j, i>=j ? vij : 0);
3127       gsl_matrix_set (lt, i, j, i<=j ? vij : 0);
3128     }
3129   }
3130
3131   /* put D into its own matrix */
3132
3133   gsl_matrix_set_zero(dm);
3134   for(i = 0; i < M; ++i) gsl_matrix_set(dm, i, i, gsl_vector_get(dv, i));
3135
3136   /* compute a = L * D * transpose(L); uses v for temp space */
3137
3138   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, dm, lt, 0.0, v);
3139   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, l,   v, 0.0, a);
3140
3141   for(i = 0; i < M; i++)
3142   {
3143     for(j = 0; j < N; j++)
3144     {
3145       const double aij = gsl_matrix_get(a, i, j);
3146       const double mij = gsl_matrix_get(m, i, j);
3147       int foo = check(aij, mij, eps);
3148       if(foo)
3149       {
3150         printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
3151       }
3152       s += foo;
3153     }
3154   }
3155
3156   gsl_vector_free(dv);
3157   gsl_matrix_free(dm);
3158   gsl_matrix_free(lt);
3159   gsl_matrix_free(l);
3160   gsl_matrix_free(v);
3161   gsl_matrix_free(a);
3162
3163   return s;
3164 }
3165
3166 int test_cholesky_decomp_unit(void)
3167 {
3168   int f;
3169   int s = 0;
3170
3171   f = test_cholesky_decomp_unit_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
3172   gsl_test(f, "  cholesky_decomp_unit hilbert(2)");
3173   s += f;
3174
3175   f = test_cholesky_decomp_unit_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
3176   gsl_test(f, "  cholesky_decomp_unit hilbert(3)");
3177   s += f;
3178
3179   f = test_cholesky_decomp_unit_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
3180   gsl_test(f, "  cholesky_decomp_unit hilbert(4)");
3181   s += f;
3182
3183   f = test_cholesky_decomp_unit_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
3184   gsl_test(f, "  cholesky_decomp_unit hilbert(12)");
3185   s += f;
3186
3187   return s;
3188 }
3189
3190 int
3191 test_choleskyc_solve_dim(const gsl_matrix_complex * m, const gsl_vector_complex * actual, double eps)
3192 {
3193   int s = 0;
3194   unsigned long i, dim = m->size1;
3195   gsl_complex z;
3196   gsl_vector_complex * rhs = gsl_vector_complex_alloc(dim);
3197   gsl_matrix_complex * u  = gsl_matrix_complex_alloc(dim,dim);
3198   gsl_vector_complex * x = gsl_vector_complex_calloc(dim);
3199   GSL_SET_IMAG(&z, 0.0);
3200   gsl_matrix_complex_memcpy(u,m);
3201   for(i=0; i<dim; i++)
3202     {
3203       GSL_SET_REAL(&z, i + 1.0);
3204       gsl_vector_complex_set(rhs, i, z);
3205     }
3206   s += gsl_linalg_complex_cholesky_decomp(u);
3207   s += gsl_linalg_complex_cholesky_solve(u, rhs, x);
3208   for(i=0; i<dim; i++) {
3209     gsl_complex y = gsl_vector_complex_get(x, i);
3210     gsl_complex a = gsl_vector_complex_get(actual, i);
3211     int foo = check(GSL_REAL(y), GSL_REAL(a), eps);
3212     int foo2 = check(GSL_IMAG(y), GSL_IMAG(a), eps);
3213     if(foo || foo2) {
3214       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, GSL_REAL(y), GSL_REAL(a));
3215       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, GSL_IMAG(y), GSL_IMAG(a));
3216     }
3217     s += foo + foo2;
3218   }
3219   gsl_vector_complex_free(x);
3220   gsl_matrix_complex_free(u);
3221   gsl_vector_complex_free(rhs);
3222
3223   return s;
3224 } /* test_choleskyc_solve_dim() */
3225
3226 int
3227 test_choleskyc_solve(void)
3228 {
3229   double data7[] = { 66,0, 0,64, 126,63, 124,-62, 61,-61, 60,60, 0,-59,
3230                      0,-64, 65,0, 62,-124, -61,-122, -60,-60, 59,-59, -58,0,
3231                      126,-63, 62,124, 308,0, 180,-240, 59,-177, 174,58, -57,-114,
3232                      124,62, -61,122, 180,240, 299,0, 174,-58, 57,171, 56,-112,
3233                      61,61, -60,60, 59,177, 174,58, 119,0, 0,112, 55,-55,
3234                      60,-60, 59,59, 174,-58, 57,-171, 0,-112, 116,0, -54,-54,
3235                      0,59, -58,0, -57,114, 56,112, 55,55, -54,54, 60,0 };
3236   double data7_sol[] = { -0.524944196428570,0.209123883928571,
3237                          1.052873883928572,0.712444196428571,
3238                          0.117568824404762,0.443191964285714,
3239                          0.412862723214286,-0.356696428571429,
3240                          0.815931919642858,-0.265820312500000,
3241                          0.777929687500000,0.119484747023810,
3242                          1.058733258928571,-0.132087053571429 };
3243   gsl_matrix_complex_view cp7 = gsl_matrix_complex_view_array(data7, 7, 7);
3244   gsl_vector_complex_view cp7_sol = gsl_vector_complex_view_array(data7_sol, 7);
3245   int f;
3246   int s = 0;
3247
3248   f = test_choleskyc_solve_dim(&cp7.matrix, &cp7_sol.vector, 1024.0 * 1024.0 * GSL_DBL_EPSILON);
3249   gsl_test(f, "  complex_cholesky_solve complex(7)");
3250   s += f;
3251
3252   return s;
3253 } /* test_choleskyc_solve() */
3254
3255 int
3256 test_choleskyc_decomp_dim(const gsl_matrix_complex * m, double eps)
3257 {
3258   int s = 0;
3259   unsigned long i,j, M = m->size1, N = m->size2;
3260
3261   gsl_matrix_complex * v  = gsl_matrix_complex_alloc(M,N);
3262   gsl_matrix_complex * a  = gsl_matrix_complex_alloc(M,N);
3263   gsl_matrix_complex * l  = gsl_matrix_complex_alloc(M,N);
3264   gsl_matrix_complex * lh  = gsl_matrix_complex_alloc(N,N);
3265
3266   gsl_matrix_complex_memcpy(v, m);
3267   gsl_matrix_complex_set_zero(l);
3268   gsl_matrix_complex_set_zero(lh);
3269
3270   s += gsl_linalg_complex_cholesky_decomp(v);
3271
3272   /* Compute L L^H */
3273
3274   for (i = 0; i < N ; i++)
3275     {
3276       for (j = 0; j <= i; j++)
3277         {
3278           gsl_complex vij = gsl_matrix_complex_get(v, i, j);
3279           gsl_matrix_complex_set (l, i, j, vij);
3280           gsl_matrix_complex_set (lh, j, i, gsl_complex_conjugate(vij));
3281         }
3282     }
3283
3284   /* compute a = l lh */
3285   gsl_blas_zgemm (CblasNoTrans,
3286                   CblasNoTrans,
3287                   GSL_COMPLEX_ONE,
3288                   l,
3289                   lh,
3290                   GSL_COMPLEX_ZERO,
3291                   a);
3292
3293   for(i=0; i<M; i++) {
3294     for(j=0; j<N; j++) {
3295       gsl_complex aij = gsl_matrix_complex_get(a, i, j);
3296       gsl_complex mij = gsl_matrix_complex_get(m, i, j);
3297       int foo_r = check(GSL_REAL(aij), GSL_REAL(mij), eps);
3298       int foo_i = check(GSL_IMAG(aij), GSL_IMAG(mij), eps);
3299       if(foo_r || foo_i) {
3300         printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, GSL_REAL(aij), GSL_REAL(mij));
3301         printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, GSL_IMAG(aij), GSL_IMAG(mij));
3302       }
3303       s += foo_r + foo_i;
3304     }
3305   }
3306
3307   gsl_matrix_complex_free(v);
3308   gsl_matrix_complex_free(a);
3309   gsl_matrix_complex_free(l);
3310   gsl_matrix_complex_free(lh);
3311
3312   return s;
3313 }
3314
3315 int
3316 test_choleskyc_decomp(void)
3317 {
3318   int f;
3319   int s = 0;
3320   double dat3[] = { 59.75,0,       49.25,172.25, 66.75,-162.75,
3321                     49.25,-172.25, 555.5,0,      -429,-333.5,
3322                     66.75,162.75,  -429,333.5,   536.5,0 };
3323   gsl_matrix_complex_view p3 = gsl_matrix_complex_view_array(dat3, 3, 3);
3324
3325   f = test_choleskyc_decomp_dim(&p3.matrix, 2 * 8.0 * GSL_DBL_EPSILON);
3326   gsl_test(f, "  complex_cholesky_decomp complex(3)");
3327   s += f;
3328
3329   return s;
3330 }
3331
3332 int
3333 test_HH_solve_dim(const gsl_matrix * m, const double * actual, double eps)
3334 {
3335   int s = 0;
3336   unsigned long i, dim = m->size1;
3337
3338   gsl_permutation * perm = gsl_permutation_alloc(dim);
3339   gsl_matrix * hh  = gsl_matrix_alloc(dim,dim);
3340   gsl_vector * d = gsl_vector_alloc(dim);
3341   gsl_vector * x = gsl_vector_alloc(dim);
3342   gsl_matrix_memcpy(hh,m);
3343   for(i=0; i<dim; i++) gsl_vector_set(x, i, i+1.0);
3344   s += gsl_linalg_HH_svx(hh, x);
3345   for(i=0; i<dim; i++) {
3346     int foo = check(gsl_vector_get(x, i),actual[i],eps);
3347     if( foo) {
3348       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
3349     }
3350     s += foo;
3351   }
3352   gsl_vector_free(x);
3353   gsl_vector_free(d);
3354   gsl_matrix_free(hh);
3355   gsl_permutation_free(perm);
3356
3357   return s;
3358 }
3359
3360 int test_HH_solve(void)
3361 {
3362   int f;
3363   int s = 0;
3364
3365   f = test_HH_solve_dim(hilb2, hilb2_solution, 8.0 * GSL_DBL_EPSILON);
3366   gsl_test(f, "  HH_solve hilbert(2)");
3367   s += f;
3368
3369   f = test_HH_solve_dim(hilb3, hilb3_solution, 128.0 * GSL_DBL_EPSILON);
3370   gsl_test(f, "  HH_solve hilbert(3)");
3371   s += f;
3372
3373   f = test_HH_solve_dim(hilb4, hilb4_solution, 2.0 * 1024.0 * GSL_DBL_EPSILON);
3374   gsl_test(f, "  HH_solve hilbert(4)");
3375   s += f;
3376
3377   f = test_HH_solve_dim(hilb12, hilb12_solution, 0.5);
3378   gsl_test(f, "  HH_solve hilbert(12)");
3379   s += f;
3380
3381   f = test_HH_solve_dim(vander2, vander2_solution, 8.0 * GSL_DBL_EPSILON);
3382   gsl_test(f, "  HH_solve vander(2)");
3383   s += f;
3384
3385   f = test_HH_solve_dim(vander3, vander3_solution, 64.0 * GSL_DBL_EPSILON);
3386   gsl_test(f, "  HH_solve vander(3)");
3387   s += f;
3388
3389   f = test_HH_solve_dim(vander4, vander4_solution, 1024.0 * GSL_DBL_EPSILON);
3390   gsl_test(f, "  HH_solve vander(4)");
3391   s += f;
3392
3393   f = test_HH_solve_dim(vander12, vander12_solution, 0.05);
3394   gsl_test(f, "  HH_solve vander(12)");
3395   s += f;
3396
3397   return s;
3398 }
3399
3400
3401 int
3402 test_TDS_solve_dim(unsigned long dim, double d, double od, const double * actual, double eps)
3403 {
3404   int s = 0;
3405   unsigned long i;
3406
3407   gsl_vector * offdiag = gsl_vector_alloc(dim-1);
3408   gsl_vector * diag = gsl_vector_alloc(dim);
3409   gsl_vector * rhs = gsl_vector_alloc(dim);
3410   gsl_vector * x = gsl_vector_alloc(dim);
3411
3412   for(i=0; i<dim; i++) {
3413     gsl_vector_set(diag, i, d);
3414     gsl_vector_set(rhs,  i, i + 1.0);
3415   }
3416   for(i=0; i<dim-1; i++) {
3417     gsl_vector_set(offdiag, i, od);
3418   }
3419
3420   s += gsl_linalg_solve_symm_tridiag(diag, offdiag, rhs, x);
3421
3422   for(i=0; i<dim; i++) {
3423     double si = gsl_vector_get(x, i);
3424     double ai = actual[i];
3425     int foo = check(si, ai, eps);
3426     if(foo) {
3427       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
3428     }
3429     s += foo;
3430   }
3431
3432   gsl_vector_free(x);
3433   gsl_vector_free(rhs);
3434   gsl_vector_free(diag);
3435   gsl_vector_free(offdiag);
3436
3437   return s;
3438 }
3439
3440
3441 int test_TDS_solve(void)
3442 {
3443   int f;
3444   int s = 0;
3445
3446   {
3447     double actual[] =  {0.0, 2.0};
3448     f = test_TDS_solve_dim(2, 1.0, 0.5, actual, 8.0 * GSL_DBL_EPSILON);
3449     gsl_test(f, "  solve_TDS dim=2 A");
3450     s += f;
3451   }
3452
3453   {
3454     double actual[] =  {3.0/8.0, 15.0/8.0};
3455     f = test_TDS_solve_dim(2, 1.0, 1.0/3.0, actual, 8.0 * GSL_DBL_EPSILON);
3456     gsl_test(f, "  solve_TDS dim=2 B");
3457     s += f;
3458   }
3459
3460   {
3461     double actual[] =  {5.0/8.0, 9.0/8.0, 2.0, 15.0/8.0, 35.0/8.0};
3462     f = test_TDS_solve_dim(5, 1.0, 1.0/3.0, actual, 8.0 * GSL_DBL_EPSILON);
3463     gsl_test(f, "  solve_TDS dim=5");
3464     s += f;
3465   }
3466
3467   return s;
3468 }
3469
3470 int
3471 test_TDS_cyc_solve_one(const unsigned long dim, const double * d, const double * od,
3472                       const double * r, const double * actual, double eps)
3473 {
3474   int s = 0;
3475   unsigned long i;
3476
3477   gsl_vector * offdiag = gsl_vector_alloc(dim);
3478   gsl_vector * diag = gsl_vector_alloc(dim);
3479   gsl_vector * rhs = gsl_vector_alloc(dim);
3480   gsl_vector * x = gsl_vector_alloc(dim);
3481
3482   for(i=0; i<dim; i++) {
3483     gsl_vector_set(diag, i, d[i]);
3484     gsl_vector_set(rhs,  i, r[i]);
3485     gsl_vector_set(offdiag, i, od[i]);
3486   }
3487
3488   s += gsl_linalg_solve_symm_cyc_tridiag(diag, offdiag, rhs, x);
3489
3490   for(i=0; i<dim; i++) {
3491     double si = gsl_vector_get(x, i);
3492     double ai = actual[i];
3493     int foo = check(si, ai, eps);
3494     if(foo) {
3495       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
3496     }
3497     s += foo;
3498   }
3499
3500   gsl_vector_free(x);
3501   gsl_vector_free(rhs);
3502   gsl_vector_free(diag);
3503   gsl_vector_free(offdiag);
3504
3505   return s;
3506 }
3507
3508 int test_TDS_cyc_solve(void)
3509 {
3510   int f;
3511   int s = 0;
3512
3513 #ifdef SUPPORT_UNDERSIZE_CYC
3514   {
3515     unsigned long dim = 1;
3516     double diag[] = {  2 };
3517     double offdiag[] = { 3 };
3518     double rhs[] = { 7 };
3519     double actual[] = { 3.5 };
3520
3521     f = test_TDS_cyc_solve_one(dim, diag, offdiag, rhs, actual, 28.0 * GSL_DBL_EPSILON);
3522     gsl_test(f, "  solve_TDS_cyc dim=%lu A", dim);
3523     s += f;
3524   }
3525
3526   {
3527     unsigned long dim = 2;
3528     double diag[] = { 1, 2 };
3529     double offdiag[] = { 3, 4 };
3530     double rhs[] = { 7, -7 };
3531     double actual[] = { -5, 4 };
3532
3533     f = test_TDS_cyc_solve_one(dim, diag, offdiag, rhs, actual, 28.0 * GSL_DBL_EPSILON);
3534     gsl_test(f, "  solve_TDS_cyc dim=%lu A", dim);
3535     s += f;
3536   }
3537 #endif
3538
3539   {
3540     unsigned long dim = 3;
3541     double diag[] = { 1, 1, 1 };
3542     double offdiag[] = { 3, 3, 3 };
3543     double rhs[] = { 7, -7, 7 };
3544     double actual[] = { -2, 5, -2 };
3545
3546     f = test_TDS_cyc_solve_one(dim, diag, offdiag, rhs, actual, 28.0 * GSL_DBL_EPSILON);
3547     gsl_test(f, "  solve_TDS_cyc dim=%lu A", dim);
3548     s += f;
3549   }
3550
3551   {
3552     unsigned long dim = 5;
3553     double diag[] = { 4, 2, 1, 2, 4 };
3554     double offdiag[] = { 1, 1, 1, 1, 1 };
3555     double rhs[] = { 30, -24, 3, 21, -30 };
3556     double actual[] = { 12, 3, -42, 42, -21 };
3557
3558     /*  f = test_TDS_cyc_solve_one(dim, diag, offdiag, rhs, actual, 7.0 * GSL_DBL_EPSILON);
3559         FIXME: bad accuracy */
3560     f = test_TDS_cyc_solve_one(dim, diag, offdiag, rhs, actual, 35.0 * GSL_DBL_EPSILON);
3561     gsl_test(f, "  solve_TDS_cyc dim=%lu B", dim);
3562     s += f;
3563   }
3564
3565   return s;
3566 }
3567
3568 int
3569 test_TDN_solve_dim(unsigned long dim, double d, double a, double b, const double * actual, double eps)
3570 {
3571   int s = 0;
3572   unsigned long i;
3573
3574   gsl_vector * abovediag = gsl_vector_alloc(dim-1);
3575   gsl_vector * belowdiag = gsl_vector_alloc(dim-1);
3576   gsl_vector * diag = gsl_vector_alloc(dim);
3577   gsl_vector * rhs = gsl_vector_alloc(dim);
3578   gsl_vector * x = gsl_vector_alloc(dim);
3579
3580   for(i=0; i<dim; i++) {
3581     gsl_vector_set(diag, i, d);
3582     gsl_vector_set(rhs,  i, i + 1.0);
3583   }
3584   for(i=0; i<dim-1; i++) {
3585     gsl_vector_set(abovediag, i, a);
3586     gsl_vector_set(belowdiag, i, b);
3587   }
3588
3589   s += gsl_linalg_solve_tridiag(diag, abovediag, belowdiag, rhs, x);
3590
3591   for(i=0; i<dim; i++) {
3592     double si = gsl_vector_get(x, i);
3593     double ai = actual[i];
3594     int foo = check(si, ai, eps);
3595     if(foo) {
3596       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
3597     }
3598     s += foo;
3599   }
3600
3601   gsl_vector_free(x);
3602   gsl_vector_free(rhs);
3603   gsl_vector_free(diag);
3604   gsl_vector_free(abovediag);
3605   gsl_vector_free(belowdiag);
3606
3607   return s;
3608 }
3609
3610
3611 int test_TDN_solve(void)
3612 {
3613   int f;
3614   int s = 0;
3615   double actual[16];
3616
3617   actual[0] =  -7.0/3.0;
3618   actual[1] =  5.0/3.0;
3619   actual[2] =  4.0/3.0;
3620   f = test_TDN_solve_dim(3, 1.0, 2.0, 1.0, actual, 2.0 * GSL_DBL_EPSILON);
3621   gsl_test(f, "  solve_TDN dim=2 A");
3622   s += f;
3623
3624   actual[0] =  0.75;
3625   actual[1] =  0.75;
3626   actual[2] =  2.625;
3627   f = test_TDN_solve_dim(3, 1.0, 1.0/3.0, 1.0/2.0, actual, 2.0 * GSL_DBL_EPSILON);
3628   gsl_test(f, "  solve_TDN dim=2 B");
3629   s += f;
3630
3631   actual[0] =  99.0/140.0;
3632   actual[1] =  41.0/35.0;
3633   actual[2] =  19.0/10.0;
3634   actual[3] =  72.0/35.0;
3635   actual[4] =  139.0/35.0;
3636   f = test_TDN_solve_dim(5, 1.0, 1.0/4.0, 1.0/2.0, actual, 35.0/8.0 * GSL_DBL_EPSILON);
3637   gsl_test(f, "  solve_TDN dim=5");
3638   s += f;
3639
3640   return s;
3641 }
3642
3643 int
3644 test_TDN_cyc_solve_dim(unsigned long dim, double d, double a, double b, const double * actual, double eps)
3645 {
3646   int s = 0;
3647   unsigned long i;
3648
3649   gsl_vector * abovediag = gsl_vector_alloc(dim);
3650   gsl_vector * belowdiag = gsl_vector_alloc(dim);
3651   gsl_vector * diag = gsl_vector_alloc(dim);
3652   gsl_vector * rhs = gsl_vector_alloc(dim);
3653   gsl_vector * x = gsl_vector_alloc(dim);
3654
3655   for(i=0; i<dim; i++) {
3656     gsl_vector_set(diag, i, d);
3657     gsl_vector_set(rhs,  i, i + 1.0);
3658   }
3659   for(i=0; i<dim; i++) {
3660     gsl_vector_set(abovediag, i, a);
3661     gsl_vector_set(belowdiag, i, b);
3662   }
3663
3664   s += gsl_linalg_solve_cyc_tridiag(diag, abovediag, belowdiag, rhs, x);
3665
3666   for(i=0; i<dim; i++) {
3667     double si = gsl_vector_get(x, i);
3668     double ai = actual[i];
3669     int foo = check(si, ai, eps);
3670     if(foo) {
3671       printf("%3lu[%lu]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
3672     }
3673     s += foo;
3674   }
3675
3676   gsl_vector_free(x);
3677   gsl_vector_free(rhs);
3678   gsl_vector_free(diag);
3679   gsl_vector_free(abovediag);
3680   gsl_vector_free(belowdiag);
3681
3682   return s;
3683 }
3684
3685
3686 int test_TDN_cyc_solve(void)
3687 {
3688   int f;
3689   int s = 0;
3690   double actual[16];
3691
3692   actual[0] =  3.0/2.0;
3693   actual[1] = -1.0/2.0;
3694   actual[2] =  1.0/2.0;
3695   f = test_TDN_cyc_solve_dim(3, 1.0, 2.0, 1.0, actual, 32.0 * GSL_DBL_EPSILON);
3696   gsl_test(f, "  solve_TDN_cyc dim=2 A");
3697   s += f;
3698
3699   actual[0] = -5.0/22.0;
3700   actual[1] = -3.0/22.0;
3701   actual[2] =  29.0/22.0;
3702   actual[3] = -9.0/22.0;
3703   actual[4] =  43.0/22.0;
3704   f = test_TDN_cyc_solve_dim(5, 3.0, 2.0, 1.0, actual, 66.0 * GSL_DBL_EPSILON);
3705   gsl_test(f, "  solve_TDN_cyc dim=5");
3706   s += f;
3707
3708   return s;
3709 }
3710
3711 int
3712 test_bidiag_decomp_dim(const gsl_matrix * m, double eps)
3713 {
3714   int s = 0;
3715   unsigned long i,j,k,r, M = m->size1, N = m->size2;
3716
3717   gsl_matrix * A  = gsl_matrix_alloc(M,N);
3718   gsl_matrix * a  = gsl_matrix_alloc(M,N);
3719   gsl_matrix * b  = gsl_matrix_alloc(N,N);
3720
3721   gsl_matrix * u  = gsl_matrix_alloc(M,N);
3722   gsl_matrix * v  = gsl_matrix_alloc(N,N);
3723
3724   gsl_vector * tau1  = gsl_vector_alloc(N);
3725   gsl_vector * tau2  = gsl_vector_alloc(N-1);
3726   gsl_vector * d  = gsl_vector_alloc(N);
3727   gsl_vector * sd  = gsl_vector_alloc(N-1);
3728
3729   gsl_matrix_memcpy(A,m);
3730
3731   s += gsl_linalg_bidiag_decomp(A, tau1, tau2);
3732   s += gsl_linalg_bidiag_unpack(A, tau1, u, tau2, v, d, sd);
3733
3734   gsl_matrix_set_zero(b);
3735   for (i = 0; i < N; i++) gsl_matrix_set(b, i,i, gsl_vector_get(d,i));
3736   for (i = 0; i < N-1; i++) gsl_matrix_set(b, i,i+1, gsl_vector_get(sd,i));
3737
3738   /* Compute A = U B V^T */
3739
3740   for (i = 0; i < M ; i++)
3741     {
3742       for (j = 0; j < N; j++)
3743         {
3744           double sum = 0;
3745
3746           for (k = 0; k < N; k++)
3747             {
3748               for (r = 0; r < N; r++)
3749                 {
3750                   sum += gsl_matrix_get(u, i, k) * gsl_matrix_get (b, k, r)
3751                     * gsl_matrix_get(v, j, r);
3752                 }
3753             }
3754           gsl_matrix_set (a, i, j, sum);
3755         }
3756     }
3757
3758   for(i=0; i<M; i++) {
3759     for(j=0; j<N; j++) {
3760       double aij = gsl_matrix_get(a, i, j);
3761       double mij = gsl_matrix_get(m, i, j);
3762       int foo = check(aij, mij, eps);
3763       if(foo) {
3764         printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
3765       }
3766       s += foo;
3767     }
3768   }
3769
3770   gsl_matrix_free(A);
3771   gsl_matrix_free(a);
3772   gsl_matrix_free(u);
3773   gsl_matrix_free(v);
3774   gsl_matrix_free(b);
3775   gsl_vector_free(tau1);
3776   gsl_vector_free(tau2);
3777   gsl_vector_free(d);
3778   gsl_vector_free(sd);
3779
3780   return s;
3781 }
3782
3783 int test_bidiag_decomp(void)
3784 {
3785   int f;
3786   int s = 0;
3787
3788   f = test_bidiag_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
3789   gsl_test(f, "  bidiag_decomp m(5,3)");
3790   s += f;
3791
3792   f = test_bidiag_decomp_dim(m97, 2 * 64.0 * GSL_DBL_EPSILON);
3793   gsl_test(f, "  bidiag_decomp m(9,7)");
3794   s += f;
3795
3796   f = test_bidiag_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
3797   gsl_test(f, "  bidiag_decomp hilbert(2)");
3798   s += f;
3799
3800   f = test_bidiag_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
3801   gsl_test(f, "  bidiag_decomp hilbert(3)");
3802   s += f;
3803
3804   f = test_bidiag_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
3805   gsl_test(f, "  bidiag_decomp hilbert(4)");
3806   s += f;
3807
3808   f = test_bidiag_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
3809   gsl_test(f, "  bidiag_decomp hilbert(12)");
3810   s += f;
3811
3812   return s;
3813 }
3814
3815 void
3816 my_error_handler (const char *reason, const char *file, int line, int err)
3817 {
3818   if (0) printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err) ;
3819 }
3820
3821 int main(void)
3822 {
3823   gsl_ieee_env_setup ();
3824   gsl_set_error_handler (&my_error_handler);
3825
3826   m11 = create_general_matrix(1,1);
3827   m51 = create_general_matrix(5,1);
3828
3829   m35 = create_general_matrix(3,5);
3830   m53 = create_general_matrix(5,3);
3831   m97 = create_general_matrix(9,7);
3832
3833   s35 = create_singular_matrix(3,5);
3834   s53 = create_singular_matrix(5,3);
3835
3836   hilb2 = create_hilbert_matrix(2);
3837   hilb3 = create_hilbert_matrix(3);
3838   hilb4 = create_hilbert_matrix(4);
3839   hilb12 = create_hilbert_matrix(12);
3840
3841   vander2 = create_vandermonde_matrix(2);
3842   vander3 = create_vandermonde_matrix(3);
3843   vander4 = create_vandermonde_matrix(4);
3844   vander12 = create_vandermonde_matrix(12);
3845
3846   moler10 = create_moler_matrix(10);
3847
3848   c7 = create_complex_matrix(7);
3849
3850   row3 = create_row_matrix(3,3);
3851   row5 = create_row_matrix(5,5);
3852   row12 = create_row_matrix(12,12);
3853
3854   A22 = create_2x2_matrix (0.0, 0.0, 0.0, 0.0);
3855   A33 = gsl_matrix_alloc(3,3);
3856   A44 = gsl_matrix_alloc(4,4);
3857   A55 = gsl_matrix_alloc(5,5);
3858
3859   inf5 = create_diagonal_matrix (inf5_data, 5);
3860   gsl_matrix_set(inf5, 3, 3, GSL_POSINF);
3861
3862   nan5 = create_diagonal_matrix (inf5_data, 5);
3863   gsl_matrix_set(nan5, 3, 3, GSL_NAN);
3864
3865   dblmin3 = create_general_matrix (3, 3);
3866   gsl_matrix_scale(dblmin3, GSL_DBL_MIN);
3867
3868   dblmin5 = create_general_matrix (5, 5);
3869   gsl_matrix_scale(dblmin5, GSL_DBL_MIN);
3870
3871   /* Matmult now obsolete */
3872 #ifdef MATMULT
3873   gsl_test(test_matmult(),               "Matrix Multiply");
3874   gsl_test(test_matmult_mod(),           "Matrix Multiply with Modification");
3875 #endif
3876   gsl_test(test_bidiag_decomp(),         "Bidiagonal Decomposition");
3877   gsl_test(test_LU_solve(),              "LU Decomposition and Solve");
3878   gsl_test(test_LUc_solve(),             "Complex LU Decomposition and Solve");
3879   gsl_test(test_QR_decomp(),             "QR Decomposition");
3880   gsl_test(test_QR_solve(),              "QR Solve");
3881   gsl_test(test_LQ_solve(),              "LQ Solve");
3882   gsl_test(test_PTLQ_solve(),            "PTLQ Solve");
3883
3884   gsl_test(test_LQ_decomp(),             "LQ Decomposition");
3885   gsl_test(test_LQ_LQsolve(),            "LQ LQ Solve");
3886   gsl_test(test_LQ_lssolve(),            "LQ LS Solve");
3887   gsl_test(test_LQ_update(),             "LQ Rank-1 Update");
3888   gsl_test(test_QRPT_decomp(),           "PTLQ Decomposition");
3889   gsl_test(test_PTLQ_solve(),            "PTLQ Solve");
3890
3891   gsl_test(test_QR_QRsolve(),            "QR QR Solve");
3892   gsl_test(test_QR_lssolve(),            "QR LS Solve");
3893   gsl_test(test_QR_update(),             "QR Rank-1 Update");
3894   gsl_test(test_QRPT_decomp(),           "QRPT Decomposition");
3895   gsl_test(test_QRPT_solve(),            "QRPT Solve");
3896   gsl_test(test_QRPT_QRsolve(),          "QRPT QR Solve");
3897   gsl_test(test_SV_decomp(),             "Singular Value Decomposition");
3898   gsl_test(test_SV_decomp_jacobi(),        "Singular Value Decomposition (Jacobi)");
3899   gsl_test(test_SV_decomp_mod(),         "Singular Value Decomposition (Mod)");
3900   gsl_test(test_SV_solve(),              "SVD Solve");
3901   gsl_test(test_cholesky_decomp(),       "Cholesky Decomposition");
3902   gsl_test(test_cholesky_decomp_unit(),  "Cholesky Decomposition [unit triangular]");
3903   gsl_test(test_cholesky_solve(),        "Cholesky Solve");
3904   gsl_test(test_choleskyc_decomp(),      "Complex Cholesky Decomposition");
3905   gsl_test(test_choleskyc_solve(),       "Complex Cholesky Solve");
3906   gsl_test(test_HH_solve(),              "Householder solve");
3907   gsl_test(test_TDS_solve(),             "Tridiagonal symmetric solve");
3908   gsl_test(test_TDS_cyc_solve(),         "Tridiagonal symmetric cyclic solve");
3909   gsl_test(test_TDN_solve(),             "Tridiagonal nonsymmetric solve");
3910   gsl_test(test_TDN_cyc_solve(),         "Tridiagonal nonsymmetric cyclic solve");
3911
3912   gsl_matrix_free(m11);
3913   gsl_matrix_free(m35);
3914   gsl_matrix_free(m51);
3915   gsl_matrix_free(m53);
3916   gsl_matrix_free(m97);
3917   gsl_matrix_free(s35);
3918   gsl_matrix_free(s53);
3919
3920   gsl_matrix_free(hilb2);
3921   gsl_matrix_free(hilb3);
3922   gsl_matrix_free(hilb4);
3923   gsl_matrix_free(hilb12);
3924
3925   gsl_matrix_free(vander2);
3926   gsl_matrix_free(vander3);
3927   gsl_matrix_free(vander4);
3928   gsl_matrix_free(vander12);
3929
3930   gsl_matrix_free(moler10);
3931
3932   gsl_matrix_complex_free(c7);
3933   gsl_matrix_free(row3);
3934   gsl_matrix_free(row5);
3935   gsl_matrix_free(row12);
3936
3937   gsl_matrix_free(A22);
3938   gsl_matrix_free(A33);
3939   gsl_matrix_free(A44);
3940   gsl_matrix_free(A55);
3941
3942   gsl_matrix_free (inf5);
3943   gsl_matrix_free (nan5);
3944
3945   exit (gsl_test_summary());
3946 }