htswanalysis/MACS/lib/gsl/gsl-1.11/matrix/test_complex_source.c

   1 /* matrix/test_complex_source.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 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 void FUNCTION (test, func) (void);
  21 void FUNCTION (test, trap) (void);
  22 void FUNCTION (test, text) (void);
  23 void FUNCTION (test, binary) (void);
  24 void FUNCTION (test, arith) (void);
  25
  26 #define TEST(expr,desc) gsl_test((expr), NAME(gsl_matrix) desc " M=%d, N=%d", M, N)
  27
  28 void
  29 FUNCTION (test, func) (void)
  30 {
  31
  32   size_t i, j;
  33   int k = 0;
  34
  35   TYPE (gsl_matrix) * m = FUNCTION (gsl_matrix, alloc) (M, N);
  36
  37   gsl_test (m->data == 0, NAME (gsl_matrix) "_alloc returns valid pointer");
  38   gsl_test (m->size1 != M, NAME (gsl_matrix) "_alloc returns valid size1");
  39   gsl_test (m->size2 != N, NAME (gsl_matrix) "_alloc returns valid size2");
  40   gsl_test (m->tda != N, NAME (gsl_matrix) "_alloc returns valid tda");
  41
  42   for (i = 0; i < M; i++)
  43     {
  44       for (j = 0; j < N; j++)
  45         {
  46           BASE z = ZERO;
  47           k++;
  48           GSL_REAL (z) = (ATOMIC) k;
  49           GSL_IMAG (z) = (ATOMIC) (k + 1000);
  50           FUNCTION (gsl_matrix, set) (m, i, j, z);
  51         }
  52     }
  53
  54   status = 0;
  55   k = 0;
  56   for (i = 0; i < M; i++)
  57     {
  58       for (j = 0; j < N; j++)
  59         {
  60           k++;
  61           if (m->data[2 * (i * N + j)] != k ||
  62               m->data[2 * (i * N + j) + 1] != k + 1000)
  63             status = 1;
  64         }
  65     }
  66
  67   gsl_test (status, NAME (gsl_matrix) "_set writes into array");
  68
  69   status = 0;
  70   k = 0;
  71   for (i = 0; i < M; i++)
  72     {
  73       for (j = 0; j < N; j++)
  74         {
  75           BASE z = FUNCTION (gsl_matrix, get) (m, i, j);
  76           k++;
  77           if (GSL_REAL (z) != k || GSL_IMAG (z) != k + 1000)
  78             status = 1;
  79         }
  80     }
  81   gsl_test (status, NAME (gsl_matrix) "_get reads from array");
  82
  83   FUNCTION (gsl_matrix, free) (m);      /* free whatever is in m */
  84
  85   m = FUNCTION (gsl_matrix, calloc) (M, N);
  86
  87   {
  88     int status = (FUNCTION(gsl_matrix,isnull)(m) != 1);
  89     TEST (status, "_isnull" DESC " on calloc matrix");
  90
  91     status = (FUNCTION(gsl_matrix,ispos)(m) != 0);
  92     TEST (status, "_ispos" DESC " on calloc matrix");
  93
  94     status = (FUNCTION(gsl_matrix,isneg)(m) != 0);
  95     TEST (status, "_isneg" DESC " on calloc matrix");
  96   }
  97
  98   for (i = 0; i < M; i++)
  99     {
 100       for (j = 0; j < N; j++)
 101         {
 102           BASE z = ZERO;
 103           FUNCTION (gsl_matrix, set) (m, i, j, z);
 104         }
 105     }
 106
 107   {
 108     status = (FUNCTION(gsl_matrix,isnull)(m) != 1);
 109     TEST (status, "_isnull" DESC " on null matrix") ;
 110
 111     status = (FUNCTION(gsl_matrix,ispos)(m) != 0);
 112     TEST (status, "_ispos" DESC " on null matrix") ;
 113
 114     status = (FUNCTION(gsl_matrix,isneg)(m) != 0);
 115     TEST (status, "_isneg" DESC " on null matrix") ;
 116   }
 117
 118
 119   k = 0;
 120   for (i = 0; i < M; i++)
 121     {
 122       for (j = 0; j < N; j++)
 123         {
 124           BASE z = ZERO;
 125           k++;
 126           GSL_REAL (z) = (ATOMIC) (k % 10);
 127           GSL_IMAG (z) = (ATOMIC) ((k + 5) % 10);
 128           FUNCTION (gsl_matrix, set) (m, i, j, z);
 129         }
 130     }
 131
 132   {
 133     status = (FUNCTION(gsl_matrix,isnull)(m) != 0);
 134     TEST (status, "_isnull" DESC " on non-negative matrix") ;
 135
 136     status = (FUNCTION(gsl_matrix,ispos)(m) != 0);
 137     TEST (status, "_ispos" DESC " on non-negative matrix") ;
 138
 139     status = (FUNCTION(gsl_matrix,isneg)(m) != 0);
 140     TEST (status, "_isneg" DESC " on non-negative matrix") ;
 141   }
 142
 143   k = 0;
 144   for (i = 0; i < M; i++)
 145     {
 146       for (j = 0; j < N; j++)
 147         {
 148           BASE z = ZERO;
 149           k++;
 150           GSL_REAL (z) = (ATOMIC) ((k % 10) - 5);
 151           GSL_IMAG (z) = (ATOMIC) (((k + 5) % 10) - 5);
 152           FUNCTION (gsl_matrix, set) (m, i, j, z);
 153         }
 154     }
 155
 156   {
 157     status = (FUNCTION(gsl_matrix,isnull)(m) != 0);
 158     TEST (status, "_isnull" DESC " on mixed matrix") ;
 159
 160     status = (FUNCTION(gsl_matrix,ispos)(m) != 0);
 161     TEST (status, "_ispos" DESC " on mixed matrix") ;
 162
 163     status = (FUNCTION(gsl_matrix,isneg)(m) != 0);
 164     TEST (status, "_isneg" DESC " on mixed matrix") ;
 165   }
 166
 167   k = 0;
 168   for (i = 0; i < M; i++)
 169     {
 170       for (j = 0; j < N; j++)
 171         {
 172           BASE z = ZERO;
 173           k++;
 174           GSL_REAL (z) = -(ATOMIC) (k % 10);
 175           GSL_IMAG (z) = -(ATOMIC) ((k + 5) % 10);
 176           FUNCTION (gsl_matrix, set) (m, i, j, z);
 177         }
 178     }
 179
 180   {
 181     status = (FUNCTION(gsl_matrix,isnull)(m) != 0);
 182     TEST (status, "_isnull" DESC " on non-positive matrix") ;
 183
 184     status = (FUNCTION(gsl_matrix,ispos)(m) != 0);
 185     TEST (status, "_ispos" DESC " on non-positive matrix") ;
 186
 187     status = (FUNCTION(gsl_matrix,isneg)(m) != 0);
 188     TEST (status, "_isneg" DESC " on non-positive matrix") ;
 189   }
 190
 191   k = 0;
 192   for (i = 0; i < M; i++)
 193     {
 194       for (j = 0; j < N; j++)
 195         {
 196           BASE z = ZERO;
 197           k++;
 198           GSL_REAL (z) = (ATOMIC) (k % 10 + 1);
 199           GSL_IMAG (z) = (ATOMIC) ((k + 5) % 10 + 1);
 200           FUNCTION (gsl_matrix, set) (m, i, j, z);
 201         }
 202     }
 203
 204   {
 205     status = (FUNCTION(gsl_matrix,isnull)(m) != 0);
 206     TEST (status, "_isnull" DESC " on positive matrix") ;
 207
 208     status = (FUNCTION(gsl_matrix,ispos)(m) != 1);
 209     TEST (status, "_ispos" DESC " on positive matrix") ;
 210
 211     status = (FUNCTION(gsl_matrix,isneg)(m) != 0);
 212     TEST (status, "_isneg" DESC " on positive matrix") ;
 213   }
 214
 215   k = 0;
 216   for (i = 0; i < M; i++)
 217     {
 218       for (j = 0; j < N; j++)
 219         {
 220           BASE z = ZERO;
 221           k++;
 222           GSL_REAL (z) = -(ATOMIC) (k % 10 + 1);
 223           GSL_IMAG (z) = -(ATOMIC) ((k + 5) % 10 + 1);
 224           FUNCTION (gsl_matrix, set) (m, i, j, z);
 225         }
 226     }
 227
 228   {
 229     status = (FUNCTION(gsl_matrix,isnull)(m) != 0);
 230     TEST (status, "_isnull" DESC " on negative matrix") ;
 231
 232     status = (FUNCTION(gsl_matrix,ispos)(m) != 0);
 233     TEST (status, "_ispos" DESC " on negative matrix") ;
 234
 235     status = (FUNCTION(gsl_matrix,isneg)(m) != 1);
 236     TEST (status, "_isneg" DESC " on negative matrix") ;
 237   }
 238
 239   FUNCTION (gsl_matrix, free) (m);      /* free whatever is in m */
 240 }
 241
 242 #if !(USES_LONGDOUBLE && !HAVE_PRINTF_LONGDOUBLE)
 243 void
 244 FUNCTION (test, text) (void)
 245 {
 246   TYPE (gsl_matrix) * m = FUNCTION (gsl_matrix, alloc) (M, N);
 247
 248   size_t i, j;
 249   int k = 0;
 250
 251   {
 252     FILE *f = fopen ("test.txt", "w");
 253     k = 0;
 254     for (i = 0; i < M; i++)
 255       {
 256         for (j = 0; j < N; j++)
 257           {
 258             BASE z;
 259             k++;
 260             GSL_REAL (z) = (ATOMIC) k;
 261             GSL_IMAG (z) = (ATOMIC) (k + 1000);
 262             FUNCTION (gsl_matrix, set) (m, i, j, z);
 263           }
 264       }
 265
 266     FUNCTION (gsl_matrix, fprintf) (f, m, OUT_FORMAT);
 267
 268     fclose (f);
 269   }
 270
 271   {
 272     FILE *f = fopen ("test.txt", "r");
 273     TYPE (gsl_matrix) * mm = FUNCTION (gsl_matrix, alloc) (M, N);
 274     status = 0;
 275
 276     FUNCTION (gsl_matrix, fscanf) (f, mm);
 277     k = 0;
 278     for (i = 0; i < M; i++)
 279       {
 280         for (j = 0; j < N; j++)
 281           {
 282             k++;
 283             if (mm->data[2 * (i * N + j)] != k
 284                 || mm->data[2 * (i * N + j) + 1] != k + 1000)
 285               status = 1;
 286           }
 287       }
 288
 289     gsl_test (status, NAME (gsl_matrix) "_fprintf and fscanf");
 290
 291     fclose (f);
 292     FUNCTION (gsl_matrix, free) (mm);
 293   }
 294
 295   FUNCTION (gsl_matrix, free) (m);
 296 }
 297 #endif
 298
 299 void
 300 FUNCTION (test, binary) (void)
 301 {
 302   TYPE (gsl_matrix) * m = FUNCTION (gsl_matrix, alloc) (M, N);
 303
 304   size_t i, j;
 305   int k = 0;
 306
 307   {
 308     FILE *f = fopen ("test.dat", "wb");
 309     k = 0;
 310     for (i = 0; i < M; i++)
 311       {
 312         for (j = 0; j < N; j++)
 313           {
 314             BASE z = ZERO;
 315             k++;
 316             GSL_REAL (z) = (ATOMIC) k;
 317             GSL_IMAG (z) = (ATOMIC) (k + 1000);
 318             FUNCTION (gsl_matrix, set) (m, i, j, z);
 319           }
 320       }
 321
 322     FUNCTION (gsl_matrix, fwrite) (f, m);
 323
 324     fclose (f);
 325   }
 326
 327   {
 328     FILE *f = fopen ("test.dat", "rb");
 329     TYPE (gsl_matrix) * mm = FUNCTION (gsl_matrix, alloc) (M, N);
 330     status = 0;
 331
 332     FUNCTION (gsl_matrix, fread) (f, mm);
 333     k = 0;
 334     for (i = 0; i < M; i++)
 335       {
 336         for (j = 0; j < N; j++)
 337           {
 338             k++;
 339             if (mm->data[2 * (i * N + j)] != k
 340                 || mm->data[2 * (i * N + j) + 1] != k + 1000)
 341               status = 1;
 342           }
 343       }
 344
 345     gsl_test (status, NAME (gsl_matrix) "_write and read");
 346
 347     fclose (f);
 348     FUNCTION (gsl_matrix, free) (mm);
 349   }
 350
 351   FUNCTION (gsl_matrix, free) (m);
 352 }
 353
 354 void
 355 FUNCTION (test, trap) (void)
 356 {
 357   TYPE (gsl_matrix) * mc = FUNCTION (gsl_matrix, alloc) (M, N);
 358   size_t i = 0, j = 0;
 359
 360   BASE z = { {(ATOMIC) 1.2, (ATOMIC) 3.4} };
 361   BASE z1;
 362
 363   status = 0;
 364   FUNCTION (gsl_matrix, set) (mc, i - 1, j, z);
 365   gsl_test (!status,
 366             NAME (gsl_matrix) "_set traps 1st index below lower bound");
 367
 368   status = 0;
 369   FUNCTION (gsl_matrix, set) (mc, i, j - 1, z);
 370   gsl_test (!status,
 371             NAME (gsl_matrix) "_set traps 2nd index below lower bound");
 372
 373   status = 0;
 374   FUNCTION (gsl_matrix, set) (mc, M + 1, 0, z);
 375   gsl_test (!status,
 376             NAME (gsl_matrix) "_set traps 1st index above upper bound");
 377
 378   status = 0;
 379   FUNCTION (gsl_matrix, set) (mc, 0, N + 1, z);
 380   gsl_test (!status,
 381             NAME (gsl_matrix) "_set traps 2nd index above upper bound");
 382
 383   status = 0;
 384   FUNCTION (gsl_matrix, set) (mc, M, 0, z);
 385   gsl_test (!status, NAME (gsl_matrix) "_set traps 1st index at upper bound");
 386
 387   status = 0;
 388   FUNCTION (gsl_matrix, set) (mc, 0, N, z);
 389   gsl_test (!status, NAME (gsl_matrix) "_set traps 2nd index at upper bound");
 390
 391   status = 0;
 392   z1 = FUNCTION (gsl_matrix, get) (mc, i - 1, 0);
 393   gsl_test (!status,
 394             NAME (gsl_matrix) "_get traps 1st index below lower bound");
 395   gsl_test (GSL_REAL (z1) != 0,
 396             NAME (gsl_matrix) "_get, zero real for 1st index below l.b.");
 397   gsl_test (GSL_IMAG (z1) != 0,
 398             NAME (gsl_matrix) "_get, zero imag for 1st index below l.b.");
 399
 400   status = 0;
 401   z1 = FUNCTION (gsl_matrix, get) (mc, 0, j - 1);
 402   gsl_test (!status,
 403             NAME (gsl_matrix) "_get traps 2nd index below lower bound");
 404   gsl_test (GSL_REAL (z1) != 0,
 405             NAME (gsl_matrix) "_get, zero real for 2nd index below l.b.");
 406   gsl_test (GSL_IMAG (z1) != 0,
 407             NAME (gsl_matrix) "_get, zero imag for 2nd index below l.b.");
 408
 409   status = 0;
 410   z1 = FUNCTION (gsl_matrix, get) (mc, M + 1, 0);
 411   gsl_test (!status,
 412             NAME (gsl_matrix) "_get traps 1st index above upper bound");
 413   gsl_test (GSL_REAL (z1) != 0,
 414             NAME (gsl_matrix) "_get, zero real for 1st index above u.b.");
 415   gsl_test (GSL_IMAG (z1) != 0,
 416             NAME (gsl_matrix) "_get, zero imag for 1st index above u.b.");
 417
 418   status = 0;
 419   z1 = FUNCTION (gsl_matrix, get) (mc, 0, N + 1);
 420   gsl_test (!status,
 421             NAME (gsl_matrix) "_get traps 2nd index above upper bound");
 422   gsl_test (GSL_REAL (z1) != 0,
 423             NAME (gsl_matrix) "_get, zero real for 2nd index above u.b.");
 424   gsl_test (GSL_IMAG (z1) != 0,
 425             NAME (gsl_matrix) "_get, zero imag for 2nd index above u.b.");
 426
 427   status = 0;
 428   z1 = FUNCTION (gsl_matrix, get) (mc, M, 0);
 429   gsl_test (!status, NAME (gsl_matrix) "_get traps 1st index at upper bound");
 430   gsl_test (GSL_REAL (z1) != 0,
 431             NAME (gsl_matrix) "_get, zero real for 1st index at u.b.");
 432   gsl_test (GSL_IMAG (z1) != 0,
 433             NAME (gsl_matrix) "_get, zero imag for 1st index at u.b.");
 434
 435   status = 0;
 436   z1 = FUNCTION (gsl_matrix, get) (mc, 0, N);
 437   gsl_test (!status, NAME (gsl_matrix) "_get traps 2nd index at upper bound");
 438   gsl_test (GSL_REAL (z1) != 0,
 439             NAME (gsl_matrix) "_get, zero real for 2nd index at u.b.");
 440   gsl_test (GSL_IMAG (z1) != 0,
 441             NAME (gsl_matrix) "_get, zero imag for 2nd index at u.b.");
 442
 443   FUNCTION (gsl_matrix, free) (mc);
 444 }
 445
 446
 447 void
 448 FUNCTION (test, arith) (void)
 449 {
 450
 451 #define P 8
 452 #define Q 12
 453 /* Must use smaller dimensions to prevent approximation of floats in float_mul_elements test*/
 454
 455   TYPE (gsl_matrix) * a = FUNCTION (gsl_matrix, alloc) (P, Q);
 456   TYPE (gsl_matrix) * b = FUNCTION (gsl_matrix, alloc) (P, Q);
 457   TYPE (gsl_matrix) * m = FUNCTION (gsl_matrix, alloc) (P, Q);
 458   size_t i, j;
 459   size_t k = 0;
 460
 461   size_t status = 0;
 462
 463   for (i = 0; i < P; i++)
 464     {
 465       for (j = 0; j < Q; j++)
 466         {
 467           BASE z, z1;
 468           GSL_REAL (z) = (ATOMIC) k;
 469           GSL_IMAG (z) = (ATOMIC) (k + 10);
 470           GSL_REAL (z1) = (ATOMIC) (k + 5);
 471           GSL_IMAG (z1) = (ATOMIC) (k + 20);
 472
 473           FUNCTION (gsl_matrix, set) (a, i, j, z);
 474           FUNCTION (gsl_matrix, set) (b, i, j, z1);
 475           k++;
 476         }
 477     }
 478
 479   {
 480     FUNCTION (gsl_matrix, memcpy) (m, a);
 481
 482     FUNCTION (gsl_matrix, add) (m, b);
 483
 484     k = 0;
 485     status = 0;
 486
 487     for (i = 0; i < P; i++)
 488       {
 489         for (j = 0; j < Q; j++)
 490           {
 491             BASE z = FUNCTION (gsl_matrix, get) (m, i, j);
 492             if (GSL_REAL (z) != (ATOMIC) (2 * k + 5) ||
 493                 GSL_IMAG (z) != (ATOMIC) (2 * k + 30))
 494               status = 1;
 495             k++;
 496           }
 497       }
 498     gsl_test (status, NAME (gsl_matrix) "_add matrix addition");
 499   }
 500
 501   {
 502     FUNCTION (gsl_matrix, memcpy) (m, a);
 503
 504     FUNCTION (gsl_matrix, sub) (m, b);
 505
 506     k = 0;
 507     status = 0;
 508
 509     for (i = 0; i < P; i++)
 510       {
 511         for (j = 0; j < Q; j++)
 512           {
 513             BASE z = FUNCTION (gsl_matrix, get) (m, i, j);
 514             if (GSL_REAL (z) != (ATOMIC) (-5)
 515                 || GSL_IMAG (z) != (ATOMIC) (-10))
 516               status = 1;
 517             k++;
 518           }
 519       }
 520     gsl_test (status, NAME (gsl_matrix) "_sub matrix subtraction");
 521   }
 522
 523   {
 524     FUNCTION (gsl_matrix, memcpy) (m, a);
 525
 526     FUNCTION (gsl_matrix, mul_elements) (m, b);
 527
 528     k = 0;
 529     status = 0;
 530
 531     for (i = 0; i < P; i++)
 532       {
 533         for (j = 0; j < Q; j++)
 534           {
 535             ATOMIC real = -(ATOMIC) (25 * k + 200);
 536             ATOMIC imag = (ATOMIC) (2 * k * k + 35 * k + 50);
 537             BASE z = FUNCTION (gsl_matrix, get) (m, i, j);
 538             if (fabs (GSL_REAL (z) - real) > 100 * BASE_EPSILON ||
 539                 fabs (GSL_IMAG (z) - imag) > 100 * BASE_EPSILON)
 540               {
 541                 status = 1;
 542 #ifdef DEBUG
 543                 printf ("%d\t%d\n", i, j);
 544                 printf (OUT_FORMAT "\n",
 545                         GSL_REAL (z) + (ATOMIC) (25 * (ATOMIC) k + 200));
 546                 printf (OUT_FORMAT "\n",
 547                         GSL_IMAG (z) - (ATOMIC) (2 * k * k + 35 * k + 50));
 548                 printf ("\n");
 549 #endif
 550               }
 551             k++;
 552           }
 553       }
 554     gsl_test (status, NAME (gsl_matrix) "_mul_elements multiplication");
 555   }
 556
 557
 558   {
 559     FUNCTION (gsl_matrix, memcpy) (m, a);
 560
 561     FUNCTION (gsl_matrix, div_elements) (m, b);
 562
 563     k = 0;
 564     status = 0;
 565
 566     for (i = 0; i < P; i++)
 567       {
 568         for (j = 0; j < Q; j++)
 569           {
 570             ATOMIC denom = (2 * k * k + 50 * k + 425);
 571             ATOMIC real = (ATOMIC) (2 * k * k + 35 * k + 200) / denom;
 572             ATOMIC imag = ((ATOMIC) (50) - (ATOMIC) (5 * k)) / denom;
 573             BASE z = FUNCTION (gsl_matrix, get) (m, i, j);
 574             if (fabs (GSL_REAL (z) - real) > 100 * BASE_EPSILON ||
 575                 fabs (GSL_IMAG (z) - imag) > 100 * BASE_EPSILON)
 576               {
 577 #ifdef DEBUG
 578                 printf (OUT_FORMAT "\t",
 579                         GSL_REAL (z) - (ATOMIC) (2 * k * k + 35 * k +
 580                                                  200) / denom);
 581                 printf (OUT_FORMAT "\n",
 582                         GSL_IMAG (z) - ((ATOMIC) (50) -
 583                                         (ATOMIC) (5 * k)) / denom);
 584 #endif
 585                 status = 1;
 586               }
 587             k++;
 588           }
 589       }
 590     gsl_test (status, NAME (gsl_matrix) "_div_elements division");
 591   }
 592
 593   FUNCTION (gsl_matrix, free) (a);
 594   FUNCTION (gsl_matrix, free) (b);
 595   FUNCTION (gsl_matrix, free) (m);
 596
 597 }