home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / ztorture.c < prev    next >
C/C++ Source or Header  |  1994-01-14  |  21KB  |  716 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. /*
  28.     This file contains a series of tests for the Meschach matrix
  29.     library, complex routines
  30. */
  31.  
  32. static char rcsid[] = "$Id: $";
  33.  
  34. #include    <stdio.h>
  35. #include    <math.h>
  36. #include     "zmatrix2.h"
  37. #include        "matlab.h"
  38.  
  39.  
  40. #define    errmesg(mesg)    printf("Error: %s error: line %d\n",mesg,__LINE__)
  41. #define notice(mesg)    printf("# Testing %s...\n",mesg);
  42.  
  43. /* extern    int    malloc_chain_check(); */
  44. /* #define MEMCHK() if ( malloc_chain_check(0) ) \
  45. { printf("Error in malloc chain: \"%s\", line %d\n", \
  46.      __FILE__, __LINE__); exit(0); } */
  47. #define    MEMCHK() 
  48.  
  49. /* cmp_perm -- returns 1 if pi1 == pi2, 0 otherwise */
  50. int    cmp_perm(pi1, pi2)
  51. PERM    *pi1, *pi2;
  52. {
  53.     int        i;
  54.  
  55.     if ( ! pi1 || ! pi2 )
  56.     error(E_NULL,"cmp_perm");
  57.     if ( pi1->size != pi2->size )
  58.     return 0;
  59.     for ( i = 0; i < pi1->size; i++ )
  60.     if ( pi1->pe[i] != pi2->pe[i] )
  61.         return 0;
  62.     return 1;
  63. }
  64.  
  65. /* px_rand -- generates sort-of random permutation */
  66. PERM    *px_rand(pi)
  67. PERM    *pi;
  68. {
  69.     int        i, j, k;
  70.  
  71.     if ( ! pi )
  72.     error(E_NULL,"px_rand");
  73.  
  74.     for ( i = 0; i < 3*pi->size; i++ )
  75.     {
  76.     j = (rand() >> 8) % pi->size;
  77.     k = (rand() >> 8) % pi->size;
  78.     px_transp(pi,j,k);
  79.     }
  80.  
  81.     return pi;
  82. }
  83.  
  84. #define    SAVE_FILE    "asx5213a.mat"
  85. #define    MATLAB_NAME    "alpha"
  86. char    name[81] = MATLAB_NAME;
  87.  
  88. void    main(argc, argv)
  89. int    argc;
  90. char    *argv[];
  91. {
  92.     ZVEC     *x = ZVNULL, *y = ZVNULL, *z = ZVNULL, *u = ZVNULL;
  93.     ZVEC    *diag = ZVNULL;
  94.     PERM    *pi1 = PNULL, *pi2 = PNULL, *pivot = PNULL;
  95.     ZMAT    *A = ZMNULL, *B = ZMNULL, *C = ZMNULL, *D = ZMNULL,
  96.     *Q = ZMNULL;
  97.     complex    ONE;
  98.     complex    z1, z2, z3;
  99.     Real    cond_est, s1, s2, s3;
  100.     int        i, seed;
  101.     FILE    *fp;
  102.     char    *cp;
  103.  
  104.  
  105.     mem_info_on(TRUE);
  106.  
  107.     setbuf(stdout,(char *)NULL);
  108.  
  109.     seed = 1111;
  110.     if ( argc > 2 )
  111.     {
  112.     printf("usage: %s [seed]\n",argv[0]);
  113.     exit(0);
  114.     }
  115.     else if ( argc == 2 )
  116.     sscanf(argv[1], "%d", &seed);
  117.  
  118.     /* set seed for rand() */
  119.     smrand(seed);
  120.  
  121.     /* print out version information */
  122.     m_version();
  123.  
  124.     printf("# Meschach Complex numbers & vectors torture test\n\n");
  125.     printf("# grep \"^Error\" the output for a listing of errors\n");
  126.     printf("# Don't panic if you see \"Error\" appearing; \n");
  127.     printf("# Also check the reported size of error\n");
  128.     printf("# This program uses randomly generated problems and therefore\n");
  129.     printf("# may occasionally produce ill-conditioned problems\n");
  130.     printf("# Therefore check the size of the error compared with MACHEPS\n");
  131.     printf("# If the error is within 1000*MACHEPS then don't worry\n");
  132.     printf("# If you get an error of size 0.1 or larger there is \n");
  133.     printf("# probably a bug in the code or the compilation procedure\n\n");
  134.     printf("# seed = %d\n",seed);
  135.  
  136.     printf("\n");
  137.  
  138.     mem_stat_mark(1);
  139.  
  140.     notice("complex arithmetic & special functions");
  141.  
  142.     ONE = zmake(1.0,0.0);
  143.     printf("# ONE = ");    z_output(ONE);
  144.     z1.re = mrand();    z1.im = mrand();
  145.     z2.re = mrand();    z2.im = mrand();
  146.     z3 = zadd(z1,z2);
  147.     if ( fabs(z1.re+z2.re-z3.re) + fabs(z1.im+z2.im-z3.im) > 10*MACHEPS )
  148.     errmesg("zadd");
  149.     z3 = zsub(z1,z2);
  150.     if ( fabs(z1.re-z2.re-z3.re) + fabs(z1.im-z2.im-z3.im) > 10*MACHEPS )
  151.     errmesg("zadd");
  152.     z3 = zmlt(z1,z2);
  153.     if ( fabs(z1.re*z2.re - z1.im*z2.im - z3.re) +
  154.      fabs(z1.im*z2.re + z1.re*z2.im - z3.im) > 10*MACHEPS )
  155.     errmesg("zmlt");
  156.     s1 = zabs(z1);
  157.     if ( fabs(s1*s1 - (z1.re*z1.re+z1.im*z1.im)) > 10*MACHEPS )
  158.     errmesg("zabs");
  159.     if ( zabs(zsub(z1,zmlt(z2,zdiv(z1,z2)))) > 10*MACHEPS ||
  160.      zabs(zsub(ONE,zdiv(z1,zmlt(z2,zdiv(z1,z2))))) > 10*MACHEPS )
  161.     errmesg("zdiv");
  162.  
  163.     z3 = zsqrt(z1);
  164.     if ( zabs(zsub(z1,zmlt(z3,z3))) > 10*MACHEPS )
  165.     errmesg("zsqrt");
  166.     if ( zabs(zsub(z1,zlog(zexp(z1)))) > 10*MACHEPS )
  167.     errmesg("zexp/zlog");
  168.     
  169.  
  170.     printf("# Check: MACHEPS = %g\n",MACHEPS);
  171.     /* allocate, initialise, copy and resize operations */
  172.     /* ZVEC */
  173.     notice("vector initialise, copy & resize");
  174.     x = zv_get(12);
  175.     y = zv_get(15);
  176.     z = zv_get(12);
  177.     zv_rand(x);
  178.     zv_rand(y);
  179.     z = zv_copy(x,z);
  180.     if ( zv_norm2(zv_sub(x,z,z)) >= MACHEPS )
  181.     errmesg("ZVEC copy");
  182.     zv_copy(x,y);
  183.     x = zv_resize(x,10);
  184.     y = zv_resize(y,10);
  185.     if ( zv_norm2(zv_sub(x,y,z)) >= MACHEPS )
  186.     errmesg("ZVEC copy/resize");
  187.     x = zv_resize(x,15);
  188.     y = zv_resize(y,15);
  189.     if ( zv_norm2(zv_sub(x,y,z)) >= MACHEPS )
  190.     errmesg("VZEC resize");
  191.  
  192.     /* ZMAT */
  193.     notice("matrix initialise, copy & resize");
  194.     A = zm_get(8,5);
  195.     B = zm_get(3,9);
  196.     C = zm_get(8,5);
  197.     zm_rand(A);
  198.     zm_rand(B);
  199.     C = zm_copy(A,C);
  200.     if ( zm_norm_inf(zm_sub(A,C,C)) >= MACHEPS )
  201.     errmesg("ZMAT copy");
  202.     zm_copy(A,B);
  203.     A = zm_resize(A,3,5);
  204.     B = zm_resize(B,3,5);
  205.     if ( zm_norm_inf(zm_sub(A,B,C)) >= MACHEPS )
  206.     errmesg("ZMAT copy/resize");
  207.     A = zm_resize(A,10,10);
  208.     B = zm_resize(B,10,10);
  209.     if ( zm_norm_inf(zm_sub(A,B,C)) >= MACHEPS )
  210.     errmesg("ZMAT resize");
  211.  
  212.     MEMCHK();
  213.  
  214.     /* PERM */
  215.     notice("permutation initialise, inverting & permuting vectors");
  216.     pi1 = px_get(15);
  217.     pi2 = px_get(12);
  218.     px_rand(pi1);
  219.     zv_rand(x);
  220.     px_zvec(pi1,x,z);
  221.     y = zv_resize(y,x->dim);
  222.     pxinv_zvec(pi1,z,y);
  223.     if ( zv_norm2(zv_sub(x,y,z)) >= MACHEPS )
  224.     errmesg("PERMute vector");
  225.  
  226.     /* testing catch() etc */
  227.     notice("error handling routines");
  228.     catch(E_NULL,
  229.       catchall(zv_add(ZVNULL,ZVNULL,ZVNULL);
  230.              errmesg("tracecatch() failure"),
  231.              printf("# tracecatch() caught error\n");
  232.              error(E_NULL,"main"));
  233.                  errmesg("catch() failure"),
  234.       printf("# catch() caught E_NULL error\n"));
  235.  
  236.     /* testing inner products and v_mltadd() etc */
  237.     notice("inner products and linear combinations");
  238.     u = zv_get(x->dim);
  239.     zv_rand(u);
  240.     zv_rand(x);
  241.     zv_resize(y,x->dim);
  242.     zv_rand(y);
  243.     zv_mltadd(y,x,zneg(zdiv(zin_prod(x,y),zin_prod(x,x))),z);
  244.     if ( zabs(zin_prod(x,z)) >= MACHEPS*x->dim )
  245.     errmesg("zv_mltadd()/zin_prod()");
  246.  
  247.     z1 = zneg(zdiv(zin_prod(x,y),zmake(zv_norm2(x)*zv_norm2(x),0.0)));
  248.     zv_mlt(z1,x,u);
  249.     zv_add(y,u,u);
  250.     if ( zv_norm2(zv_sub(u,z,u)) >= MACHEPS*x->dim )
  251.     {
  252.     errmesg("zv_mlt()/zv_norm2()");
  253.     printf("# error norm = %g\n", zv_norm2(u));
  254.     }
  255.  
  256. #ifdef ANSI_C
  257.     zv_linlist(u,x,z1,y,ONE,VNULL);
  258.     if ( zv_norm2(zv_sub(u,z,u)) >= MACHEPS*x->dim )
  259.     errmesg("zv_linlist()");
  260. #endif
  261. #ifdef VARARGS
  262.     zv_linlist(u,x,z1,y,ONE,VNULL);
  263.     if ( zv_norm2(zv_sub(u,z,u)) >= MACHEPS*x->dim )
  264.     errmesg("zv_linlist()");
  265. #endif
  266.  
  267.     MEMCHK();
  268.  
  269.     /* vector norms */
  270.     notice("vector norms");
  271.     x = zv_resize(x,12);
  272.     zv_rand(x);
  273.     for ( i = 0; i < x->dim; i++ )
  274.     if ( zabs(v_entry(x,i)) >= 0.7 )
  275.         v_set_val(x,i,ONE);
  276.         else
  277.         v_set_val(x,i,zneg(ONE));
  278.     s1 = zv_norm1(x);
  279.     s2 = zv_norm2(x);    
  280.     s3 = zv_norm_inf(x);
  281.     if ( fabs(s1 - x->dim) >= MACHEPS*x->dim ||
  282.      fabs(s2 - sqrt((double)(x->dim))) >= MACHEPS*x->dim ||
  283.      fabs(s3 - 1.0) >= MACHEPS )
  284.     errmesg("zv_norm1/2/_inf()");
  285.  
  286.     /* test matrix multiply etc */
  287.     notice("matrix multiply and invert");
  288.     A = zm_resize(A,10,10);
  289.     B = zm_resize(B,10,10);
  290.     zm_rand(A);
  291.     zm_inverse(A,B);
  292.     zm_mlt(A,B,C);
  293.     for ( i = 0; i < C->m; i++ )
  294.     m_set_val(C,i,i,zsub(m_entry(C,i,i),ONE));
  295.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  296.     errmesg("zm_inverse()/zm_mlt()");
  297.  
  298.     MEMCHK();
  299.  
  300.     /* ... and adjoints */
  301.     notice("adjoints and adjoint-multiplies");
  302.     zm_adjoint(A,A);    /* can do square matrices in situ */
  303.     zmam_mlt(A,B,C);
  304.     for ( i = 0; i < C->m; i++ )
  305.     m_set_val(C,i,i,zsub(m_entry(C,i,i),ONE));
  306.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  307.     errmesg("zm_adjoint()/zmam_mlt()");
  308.     zm_adjoint(A,A);
  309.     zm_adjoint(B,B);
  310.     zmma_mlt(A,B,C);
  311.     for ( i = 0; i < C->m; i++ )
  312.     m_set_val(C,i,i,zsub(m_entry(C,i,i),ONE));
  313.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  314.     errmesg("zm_adjoint()/zmma_mlt()");
  315.     zsm_mlt(zmake(3.71,2.753),B,B);
  316.     zmma_mlt(A,B,C);
  317.     for ( i = 0; i < C->m; i++ )
  318.     m_set_val(C,i,i,zsub(m_entry(C,i,i),zmake(3.71,-2.753)));
  319.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  320.     errmesg("szm_mlt()/zmma_mlt()");
  321.     zm_adjoint(B,B);
  322.     zsm_mlt(zdiv(ONE,zmake(3.71,-2.753)),B,B);
  323.  
  324.     MEMCHK();
  325.  
  326.     /* ... and matrix-vector multiplies */
  327.     notice("matrix-vector multiplies");
  328.     x = zv_resize(x,A->n);
  329.     y = zv_resize(y,A->m);
  330.     z = zv_resize(z,A->m);
  331.     u = zv_resize(u,A->n);
  332.     zv_rand(x);
  333.     zv_rand(y);
  334.     zmv_mlt(A,x,z);
  335.     z1 = zin_prod(y,z);
  336.     zvm_mlt(A,y,u);
  337.     z2 = zin_prod(u,x);
  338.     if ( zabs(zsub(z1,z2)) >= (MACHEPS*x->dim)*x->dim )
  339.     {
  340.     errmesg("zmv_mlt()/zvm_mlt()");
  341.     printf("# difference between inner products is %g\n",
  342.            zabs(zsub(z1,z2)));
  343.     }
  344.     zmv_mlt(B,z,u);
  345.     if ( zv_norm2(zv_sub(u,x,u)) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  346.     errmesg("zmv_mlt()/zvm_mlt()");
  347.  
  348.     MEMCHK();
  349.  
  350.     /* get/set row/col */
  351.     notice("getting and setting rows and cols");
  352.     x = zv_resize(x,A->n);
  353.     y = zv_resize(y,B->m);
  354.     x = zget_row(A,3,x);
  355.     y = zget_col(B,3,y);
  356.     if ( zabs(zsub(_zin_prod(x,y,0,Z_NOCONJ),ONE)) >=
  357.     MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  358.     errmesg("zget_row()/zget_col()");
  359.     zv_mlt(zmake(-1.0,0.0),x,x);
  360.     zv_mlt(zmake(-1.0,0.0),y,y);
  361.     zset_row(A,3,x);
  362.     zset_col(B,3,y);
  363.     zm_mlt(A,B,C);
  364.     for ( i = 0; i < C->m; i++ )
  365.     m_set_val(C,i,i,zsub(m_entry(C,i,i),ONE));
  366.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  367.     errmesg("zset_row()/zset_col()");
  368.  
  369.     MEMCHK();
  370.  
  371.     /* matrix norms */
  372.     notice("matrix norms");
  373.     A = zm_resize(A,11,15);
  374.     zm_rand(A);
  375.     s1 = zm_norm_inf(A);
  376.     B = zm_adjoint(A,B);
  377.     s2 = zm_norm1(B);
  378.     if ( fabs(s1 - s2) >= MACHEPS*A->m )
  379.     errmesg("zm_norm1()/zm_norm_inf()");
  380.     C = zmam_mlt(A,A,C);
  381.     z1.re = z1.im = 0.0;
  382.     for ( i = 0; i < C->m && i < C->n; i++ )
  383.     z1 = zadd(z1,m_entry(C,i,i));
  384.     if ( fabs(sqrt(z1.re) - zm_norm_frob(A)) >= MACHEPS*A->m*A->n )
  385.     errmesg("zm_norm_frob");
  386.  
  387.     MEMCHK();
  388.     
  389.     /* permuting rows and columns */
  390.     /******************************
  391.     notice("permuting rows & cols");
  392.     A = zm_resize(A,11,15);
  393.     B = zm_resize(B,11,15);
  394.     pi1 = px_resize(pi1,A->m);
  395.     px_rand(pi1);
  396.     x = zv_resize(x,A->n);
  397.     y = zmv_mlt(A,x,y);
  398.     px_rows(pi1,A,B);
  399.     px_zvec(pi1,y,z);
  400.     zmv_mlt(B,x,u);
  401.     if ( zv_norm2(zv_sub(z,u,u)) >= MACHEPS*A->m )
  402.     errmesg("px_rows()");
  403.     pi1 = px_resize(pi1,A->n);
  404.     px_rand(pi1);
  405.     px_cols(pi1,A,B);
  406.     pxinv_zvec(pi1,x,z);
  407.     zmv_mlt(B,z,u);
  408.     if ( zv_norm2(zv_sub(y,u,u)) >= MACHEPS*A->n )
  409.     errmesg("px_cols()");
  410.     ******************************/
  411.  
  412.     MEMCHK();
  413.  
  414.     /* MATLAB save/load */
  415.     notice("MATLAB save/load");
  416.     A = zm_resize(A,12,11);
  417.     if ( (fp=fopen(SAVE_FILE,"w")) == (FILE *)NULL )
  418.     printf("Cannot perform MATLAB save/load test\n");
  419.     else
  420.     {
  421.     zm_rand(A);
  422.     zm_save(fp, A, name);
  423.     fclose(fp);
  424.     if ( (fp=fopen(SAVE_FILE,"r")) == (FILE *)NULL )
  425.         printf("Cannot open save file \"%s\"\n",SAVE_FILE);
  426.     else
  427.     {
  428.         ZM_FREE(B);
  429.         B = zm_load(fp,&cp);
  430.         if ( strcmp(name,cp) || zm_norm1(zm_sub(A,B,C)) >=
  431.          MACHEPS*A->m )
  432.         {
  433.         errmesg("zm_load()/zm_save()");
  434.         printf("# orig. name = %s, restored name = %s\n", name, cp);
  435.         printf("# orig. A =\n");    zm_output(A);
  436.         printf("# restored A =\n");    zm_output(B);
  437.         }
  438.     }
  439.     }
  440.  
  441.     MEMCHK();
  442.  
  443.     /* Now, onto matrix factorisations */
  444.     A = zm_resize(A,10,10);
  445.     B = zm_resize(B,A->m,A->n);
  446.     zm_copy(A,B);
  447.     x = zv_resize(x,A->n);
  448.     y = zv_resize(y,A->m);
  449.     z = zv_resize(z,A->n);
  450.     u = zv_resize(u,A->m);
  451.     zv_rand(x);
  452.     zmv_mlt(B,x,y);
  453.     z = zv_copy(x,z);
  454.  
  455.     notice("LU factor/solve");
  456.     pivot = px_get(A->m);
  457.     zLUfactor(A,pivot);
  458.     tracecatch(zLUsolve(A,pivot,y,x),"main");
  459.     tracecatch(cond_est = zLUcondest(A,pivot),"main");
  460.     printf("# cond(A) approx= %g\n", cond_est);
  461.     if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est)
  462.     {
  463.     errmesg("zLUfactor()/zLUsolve()");
  464.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  465.            zv_norm2(zv_sub(x,z,u)), MACHEPS);
  466.     }
  467.  
  468.  
  469.     zv_copy(y,x);
  470.     tracecatch(zLUsolve(A,pivot,x,x),"main");
  471.     tracecatch(cond_est = zLUcondest(A,pivot),"main");
  472.     if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est)
  473.     {
  474.     errmesg("zLUfactor()/zLUsolve()");
  475.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  476.            zv_norm2(zv_sub(x,z,u)), MACHEPS);
  477.     }
  478.  
  479.     zvm_mlt(B,z,y);
  480.     zv_copy(y,x);
  481.     tracecatch(zLUAsolve(A,pivot,x,x),"main");
  482.     if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est)
  483.     {
  484.     errmesg("zLUfactor()/zLUAsolve()");
  485.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  486.            zv_norm2(zv_sub(x,z,u)), MACHEPS);
  487.     }
  488.  
  489.     MEMCHK();
  490.  
  491.     /* QR factorisation */
  492.     zm_copy(B,A);
  493.     zmv_mlt(B,z,y);
  494.     notice("QR factor/solve:");
  495.     diag = zv_get(A->m);
  496.     zQRfactor(A,diag);
  497.     zQRsolve(A,diag,y,x);
  498.     if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est )
  499.     {
  500.     errmesg("zQRfactor()/zQRsolve()");
  501.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  502.            zv_norm2(zv_sub(x,z,u)), MACHEPS);
  503.     }
  504.     printf("# QR cond(A) approx= %g\n", zQRcondest(A));
  505.     Q = zm_get(A->m,A->m);
  506.     zmakeQ(A,diag,Q);
  507.     zmakeR(A,A);
  508.     zm_mlt(Q,A,C);
  509.     zm_sub(B,C,C);
  510.     if ( zm_norm1(C) >= MACHEPS*zm_norm1(Q)*zm_norm1(B) )
  511.     {
  512.     errmesg("zQRfactor()/zmakeQ()/zmakeR()");
  513.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  514.            zm_norm1(C), MACHEPS);
  515.     }
  516.  
  517.     MEMCHK();
  518.  
  519.     /* now try with a non-square matrix */
  520.     A = zm_resize(A,15,7);
  521.     zm_rand(A);
  522.     B = zm_copy(A,B);
  523.     diag = zv_resize(diag,A->n);
  524.     x = zv_resize(x,A->n);
  525.     y = zv_resize(y,A->m);
  526.     zv_rand(y);
  527.     zQRfactor(A,diag);
  528.     x = zQRsolve(A,diag,y,x);
  529.     /* z is the residual vector */
  530.     zmv_mlt(B,x,z);    zv_sub(z,y,z);
  531.     /* check B*.z = 0 */
  532.     zvm_mlt(B,z,u);
  533.     if ( zv_norm2(u) >= 100*MACHEPS*zm_norm1(B)*zv_norm2(y) )
  534.     {
  535.     errmesg("zQRfactor()/zQRsolve()");
  536.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  537.            zv_norm2(u), MACHEPS);
  538.     }
  539.     Q = zm_resize(Q,A->m,A->m);
  540.     zmakeQ(A,diag,Q);
  541.     zmakeR(A,A);
  542.     zm_mlt(Q,A,C);
  543.     zm_sub(B,C,C);
  544.     if ( zm_norm1(C) >= MACHEPS*zm_norm1(Q)*zm_norm1(B) )
  545.     {
  546.     errmesg("zQRfactor()/zmakeQ()/zmakeR()");
  547.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  548.            zm_norm1(C), MACHEPS);
  549.     }
  550.     D = zm_get(A->m,Q->m);
  551.     zmam_mlt(Q,Q,D);
  552.     for ( i = 0; i < D->m; i++ )
  553.     m_set_val(D,i,i,zsub(m_entry(D,i,i),ONE));
  554.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q) )
  555.     {
  556.     errmesg("QRfactor()/makeQ()/makeR()");
  557.     printf("# QR orthogonality error = %g [cf MACHEPS = %g]\n",
  558.            zm_norm1(D), MACHEPS);
  559.     }
  560.  
  561.     MEMCHK();
  562.  
  563.     /* QRCP factorisation */
  564.     zm_copy(B,A);
  565.     notice("QR factor/solve with column pivoting");
  566.     pivot = px_resize(pivot,A->n);
  567.     zQRCPfactor(A,diag,pivot);
  568.     z = zv_resize(z,A->n);
  569.     zQRCPsolve(A,diag,pivot,y,z);
  570.     /* pxinv_zvec(pivot,z,x); */
  571.     /* now compute residual (z) vector */
  572.     zmv_mlt(B,x,z);    zv_sub(z,y,z);
  573.     /* check B^T.z = 0 */
  574.     zvm_mlt(B,z,u);
  575.     if ( zv_norm2(u) >= MACHEPS*zm_norm1(B)*zv_norm2(y) )
  576.     {
  577.     errmesg("QRCPfactor()/QRsolve()");
  578.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  579.            zv_norm2(u), MACHEPS);
  580.     }
  581.  
  582.     Q = zm_resize(Q,A->m,A->m);
  583.     zmakeQ(A,diag,Q);
  584.     zmakeR(A,A);
  585.     zm_mlt(Q,A,C);
  586.     ZM_FREE(D);
  587.     D = zm_get(B->m,B->n);
  588.     /******************************
  589.     px_cols(pivot,C,D);
  590.     zm_sub(B,D,D);
  591.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm1(B) )
  592.     {
  593.     errmesg("QRCPfactor()/makeQ()/makeR()");
  594.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  595.            zm_norm1(D), MACHEPS);
  596.     }
  597.     ******************************/
  598.  
  599.     /* Now check eigenvalue/SVD routines */
  600.     notice("complex Schur routines");
  601.     A = zm_resize(A,11,11);
  602.     B = zm_resize(B,A->m,A->n);
  603.     C = zm_resize(C,A->m,A->n);
  604.     D = zm_resize(D,A->m,A->n);
  605.     Q = zm_resize(Q,A->m,A->n);
  606.  
  607.     MEMCHK();
  608.  
  609.     /* now test complex Schur decomposition */
  610.     /* zm_copy(A,B); */
  611.     ZM_FREE(A);
  612.     A = zm_get(11,11);
  613.     zm_rand(A);
  614.     B = zm_copy(A,B);
  615.     MEMCHK();
  616.  
  617.     B = zschur(B,Q);
  618.     MEMCHK();
  619.  
  620.     zm_mlt(Q,B,C);
  621.     zmma_mlt(C,Q,D);
  622.     MEMCHK();
  623.     zm_sub(A,D,D);
  624.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q)*zm_norm1(B)*5 )
  625.     {
  626.     errmesg("zschur()");
  627.     printf("# Schur reconstruction error = %g [cf MACHEPS = %g]\n",
  628.            zm_norm1(D), MACHEPS);
  629.     }
  630.  
  631.     /* orthogonality check */
  632.     zmma_mlt(Q,Q,D);
  633.     for ( i = 0; i < D->m; i++ )
  634.     m_set_val(D,i,i,zsub(m_entry(D,i,i),ONE));
  635.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q)*10 )
  636.     {
  637.     errmesg("zschur()");
  638.     printf("# Schur orthogonality error = %g [cf MACHEPS = %g]\n",
  639.            zm_norm1(D), MACHEPS);
  640.     }
  641.  
  642.     MEMCHK();
  643.  
  644.     /* now test SVD */
  645.     /******************************
  646.     A = zm_resize(A,11,7);
  647.     zm_rand(A);
  648.     U = zm_get(A->n,A->n);
  649.     Q = zm_resize(Q,A->m,A->m);
  650.     u = zv_resize(u,max(A->m,A->n));
  651.     svd(A,Q,U,u);
  652.     ******************************/
  653.     /* check reconstruction of A */
  654.     /******************************
  655.     D = zm_resize(D,A->m,A->n);
  656.     C = zm_resize(C,A->m,A->n);
  657.     zm_zero(D);
  658.     for ( i = 0; i < min(A->m,A->n); i++ )
  659.     zm_set_val(D,i,i,v_entry(u,i));
  660.     zmam_mlt(Q,D,C);
  661.     zm_mlt(C,U,D);
  662.     zm_sub(A,D,D);
  663.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(U)*zm_norm_inf(Q)*zm_norm1(A) )
  664.     {
  665.     errmesg("svd()");
  666.     printf("# SVD reconstruction error = %g [cf MACHEPS = %g]\n",
  667.            zm_norm1(D), MACHEPS);
  668.     }
  669.     ******************************/
  670.     /* check orthogonality of Q and U */
  671.     /******************************
  672.     D = zm_resize(D,Q->n,Q->n);
  673.     zmam_mlt(Q,Q,D);
  674.     for ( i = 0; i < D->m; i++ )
  675.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  676.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q)*5 )
  677.     {
  678.     errmesg("svd()");
  679.     printf("# SVD orthognality error (Q) = %g [cf MACHEPS = %g\n",
  680.            zm_norm1(D), MACHEPS);
  681.     }
  682.     D = zm_resize(D,U->n,U->n);
  683.     zmam_mlt(U,U,D);
  684.     for ( i = 0; i < D->m; i++ )
  685.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  686.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(U)*zm_norm_inf(U)*5 )
  687.     {
  688.     errmesg("svd()");
  689.     printf("# SVD orthognality error (U) = %g [cf MACHEPS = %g\n",
  690.            zm_norm1(D), MACHEPS);
  691.     }
  692.     for ( i = 0; i < u->dim; i++ )
  693.     if ( v_entry(u,i) < 0 || (i < u->dim-1 &&
  694.                   v_entry(u,i+1) > v_entry(u,i)) )
  695.         break;
  696.     if ( i < u->dim )
  697.     {
  698.     errmesg("svd()");
  699.     printf("# SVD sorting error\n");
  700.     }
  701.     ******************************/
  702.  
  703.     ZV_FREE(x);    ZV_FREE(y);    ZV_FREE(z);
  704.     ZV_FREE(u);    ZV_FREE(diag);
  705.     PX_FREE(pi1);    PX_FREE(pi2);    PX_FREE(pivot);
  706.     ZM_FREE(A);    ZM_FREE(B);    ZM_FREE(C);
  707.     ZM_FREE(D);    ZM_FREE(Q);
  708.  
  709.     mem_stat_free(1);
  710.  
  711.     MEMCHK();
  712.     printf("# Finished torture test for complex numbers/vectors/matrices\n");
  713.     mem_info();
  714. }
  715.  
  716.