home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / torture.c < prev    next >
C/C++ Source or Header  |  1994-01-14  |  29KB  |  1,042 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.     This file contains a series of tests for the Meschach matrix
  28.     library, parts 1 and 2
  29. */
  30.  
  31. static char rcsid[] = "$Id: $";
  32.  
  33. #include    <stdio.h>
  34. #include    <math.h>
  35. #include    "matrix2.h"
  36. #include        "matlab.h"
  37.  
  38. #define    errmesg(mesg)    printf("Error: %s error: line %d\n",mesg,__LINE__)
  39. #define notice(mesg)    printf("# Testing %s...\n",mesg);
  40.  
  41. static char *test_err_list[] = {
  42.    "unknown error",            /* 0 */
  43.    "testing error messages",        /* 1 */
  44.    "unexpected end-of-file"        /* 2 */
  45. };
  46.  
  47.  
  48. #define MAX_TEST_ERR   (sizeof(test_err_list)/sizeof(char *))
  49.  
  50. /* extern    int    malloc_chain_check(); */
  51. /* #define MEMCHK() if ( malloc_chain_check(0) ) \
  52. { printf("Error in malloc chain: \"%s\", line %d\n", \
  53.      __FILE__, __LINE__); exit(0); } */
  54. #define    MEMCHK() 
  55.  
  56. /* cmp_perm -- returns 1 if pi1 == pi2, 0 otherwise */
  57. int    cmp_perm(pi1, pi2)
  58. PERM    *pi1, *pi2;
  59. {
  60.     int        i;
  61.  
  62.     if ( ! pi1 || ! pi2 )
  63.     error(E_NULL,"cmp_perm");
  64.     if ( pi1->size != pi2->size )
  65.     return 0;
  66.     for ( i = 0; i < pi1->size; i++ )
  67.     if ( pi1->pe[i] != pi2->pe[i] )
  68.         return 0;
  69.     return 1;
  70. }
  71.  
  72. /* px_rand -- generates sort-of random permutation */
  73. PERM    *px_rand(pi)
  74. PERM    *pi;
  75. {
  76.     int        i, j, k;
  77.  
  78.     if ( ! pi )
  79.     error(E_NULL,"px_rand");
  80.  
  81.     for ( i = 0; i < 3*pi->size; i++ )
  82.     {
  83.     j = (rand() >> 8) % pi->size;
  84.     k = (rand() >> 8) % pi->size;
  85.     px_transp(pi,j,k);
  86.     }
  87.  
  88.     return pi;
  89. }
  90.  
  91. #define    SAVE_FILE    "asx5213a.mat"
  92. #define    MATLAB_NAME    "alpha"
  93. char    name[81] = MATLAB_NAME;
  94.  
  95. int main(argc, argv)
  96. int    argc;
  97. char    *argv[];
  98. {
  99.    VEC    *x = VNULL, *y = VNULL, *z = VNULL, *u = VNULL, *v = VNULL, 
  100.         *w = VNULL;
  101.    VEC    *diag = VNULL, *beta = VNULL;
  102.    PERM    *pi1 = PNULL, *pi2 = PNULL, *pi3 = PNULL, *pivot = PNULL, 
  103.         *blocks = PNULL;
  104.    MAT    *A = MNULL, *B = MNULL, *C = MNULL, *D = MNULL, *Q = MNULL, 
  105.         *U = MNULL;
  106.    BAND *bA, *bB, *bC;
  107.    Real    cond_est, s1, s2, s3;
  108.    int    i, j, seed;
  109.    FILE    *fp;
  110.    char    *cp;
  111.  
  112.  
  113.     mem_info_on(TRUE);
  114.  
  115.     setbuf(stdout,(char *)NULL);
  116.  
  117.     seed = 1111;
  118.     if ( argc > 2 )
  119.     {
  120.     printf("usage: %s [seed]\n",argv[0]);
  121.     exit(0);
  122.     }
  123.     else if ( argc == 2 )
  124.     sscanf(argv[1], "%d", &seed);
  125.  
  126.     /* set seed for rand() */
  127.     smrand(seed);
  128.  
  129.     mem_stat_mark(1);
  130.  
  131.     /* print version information */
  132.     m_version();
  133.  
  134.     printf("# grep \"^Error\" the output for a listing of errors\n");
  135.     printf("# Don't panic if you see \"Error\" appearing; \n");
  136.     printf("# Also check the reported size of error\n");
  137.     printf("# This program uses randomly generated problems and therefore\n");
  138.     printf("# may occasionally produce ill-conditioned problems\n");
  139.     printf("# Therefore check the size of the error compared with MACHEPS\n");
  140.     printf("# If the error is within 1000*MACHEPS then don't worry\n");
  141.     printf("# If you get an error of size 0.1 or larger there is \n");
  142.     printf("# probably a bug in the code or the compilation procedure\n\n");
  143.     printf("# seed = %d\n",seed);
  144.  
  145.     printf("# Check: MACHEPS = %g\n",MACHEPS);
  146.     /* allocate, initialise, copy and resize operations */
  147.     /* VEC */
  148.     notice("vector initialise, copy & resize");
  149.     x = v_get(12);
  150.     y = v_get(15);
  151.     z = v_get(12);
  152.     v_rand(x);
  153.     v_rand(y);
  154.     z = v_copy(x,z);
  155.     if ( v_norm2(v_sub(x,z,z)) >= MACHEPS )
  156.     errmesg("VEC copy");
  157.     v_copy(x,y);
  158.     x = v_resize(x,10);
  159.     y = v_resize(y,10);
  160.     if ( v_norm2(v_sub(x,y,z)) >= MACHEPS )
  161.     errmesg("VEC copy/resize");
  162.     x = v_resize(x,15);
  163.     y = v_resize(y,15);
  164.     if ( v_norm2(v_sub(x,y,z)) >= MACHEPS )
  165.     errmesg("VEC resize");
  166.  
  167.     /* MAT */
  168.     notice("matrix initialise, copy & resize");
  169.     A = m_get(8,5);
  170.     B = m_get(3,9);
  171.     C = m_get(8,5);
  172.     m_rand(A);
  173.     m_rand(B);
  174.     C = m_copy(A,C);
  175.     if ( m_norm_inf(m_sub(A,C,C)) >= MACHEPS )
  176.     errmesg("MAT copy");
  177.     m_copy(A,B);
  178.     A = m_resize(A,3,5);
  179.     B = m_resize(B,3,5);
  180.     if ( m_norm_inf(m_sub(A,B,C)) >= MACHEPS )
  181.     errmesg("MAT copy/resize");
  182.     A = m_resize(A,10,10);
  183.     B = m_resize(B,10,10);
  184.     if ( m_norm_inf(m_sub(A,B,C)) >= MACHEPS )
  185.     errmesg("MAT resize");
  186.  
  187.     MEMCHK();
  188.  
  189.     /* PERM */
  190.     notice("permutation initialise, inverting & permuting vectors");
  191.     pi1 = px_get(15);
  192.     pi2 = px_get(12);
  193.     px_rand(pi1);
  194.     v_rand(x);
  195.     px_vec(pi1,x,z);
  196.     y = v_resize(y,x->dim);
  197.     pxinv_vec(pi1,z,y);
  198.     if ( v_norm2(v_sub(x,y,z)) >= MACHEPS )
  199.     errmesg("PERMute vector");
  200.     pi2 = px_inv(pi1,pi2);
  201.     pi3 = px_mlt(pi1,pi2,PNULL);
  202.     for ( i = 0; i < pi3->size; i++ )
  203.     if ( pi3->pe[i] != i )
  204.         errmesg("PERM inverse/multiply");
  205.  
  206.     /* testing catch() etc */
  207.     notice("error handling routines");
  208.     catch(E_NULL,
  209.       catchall(v_add(VNULL,VNULL,VNULL);
  210.              errmesg("tracecatch() failure"),
  211.              printf("# tracecatch() caught error\n");
  212.              error(E_NULL,"main"));
  213.                  errmesg("catch() failure"),
  214.       printf("# catch() caught E_NULL error\n"));
  215.  
  216.     /* testing attaching a new error list (error list 2) */
  217.  
  218.     notice("attaching error lists");
  219.     printf("# IT IS NOT A REAL WARNING ... \n");
  220.     err_list_attach(2,MAX_TEST_ERR,test_err_list,TRUE);
  221.     if (!err_is_list_attached(2)) 
  222.        errmesg("attaching the error list 2");
  223.     ev_err(__FILE__,1,__LINE__,"main",2);
  224.     err_list_free(2);
  225.     if (err_is_list_attached(2)) 
  226.        errmesg("detaching the error list 2");
  227.  
  228.     /* testing inner products and v_mltadd() etc */
  229.     notice("inner products and linear combinations");
  230.     u = v_get(x->dim);
  231.     v_rand(u);
  232.     v_rand(x);
  233.     v_resize(y,x->dim);
  234.     v_rand(y);
  235.     v_mltadd(y,x,-in_prod(x,y)/in_prod(x,x),z);
  236.     if ( fabs(in_prod(x,z)) >= MACHEPS*x->dim )
  237.     errmesg("v_mltadd()/in_prod()");
  238.     s1 = -in_prod(x,y)/(v_norm2(x)*v_norm2(x));
  239.     sv_mlt(s1,x,u);
  240.     v_add(y,u,u);
  241.     if ( v_norm2(v_sub(u,z,u)) >= MACHEPS*x->dim )
  242.     errmesg("sv_mlt()/v_norm2()");
  243.  
  244. #ifdef ANSI_C 
  245.     v_linlist(u,x,s1,y,1.0,VNULL);
  246.     if ( v_norm2(v_sub(u,z,u)) >= MACHEPS*x->dim )
  247.     errmesg("v_linlist()");
  248. #endif
  249. #ifdef VARARGS
  250.     v_linlist(u,x,s1,y,1.0,VNULL);
  251.     if ( v_norm2(v_sub(u,z,u)) >= MACHEPS*x->dim )
  252.     errmesg("v_linlist()");
  253. #endif
  254.  
  255.  
  256.     MEMCHK();
  257.  
  258.     /* vector norms */
  259.     notice("vector norms");
  260.     x = v_resize(x,12);
  261.     v_rand(x);
  262.     for ( i = 0; i < x->dim; i++ )
  263.     if ( v_entry(x,i) >= 0.5 )
  264.         v_set_val(x,i,1.0);
  265.         else
  266.         v_set_val(x,i,-1.0);
  267.     s1 = v_norm1(x);
  268.     s2 = v_norm2(x);    
  269.     s3 = v_norm_inf(x);
  270.     if ( fabs(s1 - x->dim) >= MACHEPS*x->dim ||
  271.      fabs(s2 - sqrt((Real)(x->dim))) >= MACHEPS*x->dim ||
  272.      fabs(s3 - 1.0) >= MACHEPS )
  273.     errmesg("v_norm1/2/_inf()");
  274.  
  275.     /* test matrix multiply etc */
  276.     notice("matrix multiply and invert");
  277.     A = m_resize(A,10,10);
  278.     B = m_resize(B,10,10);
  279.     m_rand(A);
  280.     m_inverse(A,B);
  281.     m_mlt(A,B,C);
  282.     for ( i = 0; i < C->m; i++ )
  283.     m_set_val(C,i,i,m_entry(C,i,i)-1.0);
  284.     if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  285.     errmesg("m_inverse()/m_mlt()");
  286.  
  287.     MEMCHK();
  288.  
  289.     /* ... and transposes */
  290.     notice("transposes and transpose-multiplies");
  291.     m_transp(A,A);    /* can do square matrices in situ */
  292.     mtrm_mlt(A,B,C);
  293.     for ( i = 0; i < C->m; i++ )
  294.     m_set_val(C,i,i,m_entry(C,i,i)-1.0);
  295.     if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  296.     errmesg("m_transp()/mtrm_mlt()");
  297.     m_transp(A,A);
  298.     m_transp(B,B);
  299.     mmtr_mlt(A,B,C);
  300.     for ( i = 0; i < C->m; i++ )
  301.     m_set_val(C,i,i,m_entry(C,i,i)-1.0);
  302.     if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  303.     errmesg("m_transp()/mmtr_mlt()");
  304.     sm_mlt(3.71,B,B);
  305.     mmtr_mlt(A,B,C);
  306.     for ( i = 0; i < C->m; i++ )
  307.     m_set_val(C,i,i,m_entry(C,i,i)-3.71);
  308.     if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  309.     errmesg("sm_mlt()/mmtr_mlt()");
  310.     m_transp(B,B);
  311.     sm_mlt(1.0/3.71,B,B);
  312.  
  313.     MEMCHK();
  314.  
  315.     /* ... and matrix-vector multiplies */
  316.     notice("matrix-vector multiplies");
  317.     x = v_resize(x,A->n);
  318.     y = v_resize(y,A->m);
  319.     z = v_resize(z,A->m);
  320.     u = v_resize(u,A->n);
  321.     v_rand(x);
  322.     v_rand(y);
  323.     mv_mlt(A,x,z);
  324.     s1 = in_prod(y,z);
  325.     vm_mlt(A,y,u);
  326.     s2 = in_prod(u,x);
  327.     if ( fabs(s1 - s2) >= (MACHEPS*x->dim)*x->dim )
  328.     errmesg("mv_mlt()/vm_mlt()");
  329.     mv_mlt(B,z,u);
  330.     if ( v_norm2(v_sub(u,x,u)) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  331.     errmesg("mv_mlt()/m_inverse()");
  332.  
  333.     MEMCHK();
  334.  
  335.     /* get/set row/col */
  336.     notice("getting and setting rows and cols");
  337.     x = v_resize(x,A->n);
  338.     y = v_resize(y,B->m);
  339.     x = get_row(A,3,x);
  340.     y = get_col(B,3,y);
  341.     if ( fabs(in_prod(x,y) - 1.0) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  342.     errmesg("get_row()/get_col()");
  343.     sv_mlt(-1.0,x,x);
  344.     sv_mlt(-1.0,y,y);
  345.     set_row(A,3,x);
  346.     set_col(B,3,y);
  347.     m_mlt(A,B,C);
  348.     for ( i = 0; i < C->m; i++ )
  349.     m_set_val(C,i,i,m_entry(C,i,i)-1.0);
  350.     if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  351.     errmesg("set_row()/set_col()");
  352.  
  353.     MEMCHK();
  354.  
  355.     /* matrix norms */
  356.     notice("matrix norms");
  357.     A = m_resize(A,11,15);
  358.     m_rand(A);
  359.     s1 = m_norm_inf(A);
  360.     B = m_transp(A,B);
  361.     s2 = m_norm1(B);
  362.     if ( fabs(s1 - s2) >= MACHEPS*A->m )
  363.     errmesg("m_norm1()/m_norm_inf()");
  364.     C = mtrm_mlt(A,A,C);
  365.     s1 = 0.0;
  366.     for ( i = 0; i < C->m && i < C->n; i++ )
  367.     s1 += m_entry(C,i,i);
  368.     if ( fabs(sqrt(s1) - m_norm_frob(A)) >= MACHEPS*A->m*A->n )
  369.     errmesg("m_norm_frob");
  370.  
  371.     MEMCHK();
  372.     
  373.     /* permuting rows and columns */
  374.     notice("permuting rows & cols");
  375.     A = m_resize(A,11,15);
  376.     B = m_resize(B,11,15);
  377.     pi1 = px_resize(pi1,A->m);
  378.     px_rand(pi1);
  379.     x = v_resize(x,A->n);
  380.     y = mv_mlt(A,x,y);
  381.     px_rows(pi1,A,B);
  382.     px_vec(pi1,y,z);
  383.     mv_mlt(B,x,u);
  384.     if ( v_norm2(v_sub(z,u,u)) >= MACHEPS*A->m )
  385.     errmesg("px_rows()");
  386.     pi1 = px_resize(pi1,A->n);
  387.     px_rand(pi1);
  388.     px_cols(pi1,A,B);
  389.     pxinv_vec(pi1,x,z);
  390.     mv_mlt(B,z,u);
  391.     if ( v_norm2(v_sub(y,u,u)) >= MACHEPS*A->n )
  392.     errmesg("px_cols()");
  393.  
  394.     MEMCHK();
  395.  
  396.     /* MATLAB save/load */
  397.     notice("MATLAB save/load");
  398.     A = m_resize(A,12,11);
  399.     if ( (fp=fopen(SAVE_FILE,"w")) == (FILE *)NULL )
  400.     printf("Cannot perform MATLAB save/load test\n");
  401.     else
  402.     {
  403.     m_rand(A);
  404.     m_save(fp, A, name);
  405.     fclose(fp);
  406.     if ( (fp=fopen(SAVE_FILE,"r")) == (FILE *)NULL )
  407.         printf("Cannot open save file \"%s\"\n",SAVE_FILE);
  408.     else
  409.     {
  410.         M_FREE(B);
  411.         B = m_load(fp,&cp);
  412.         if ( strcmp(name,cp) || m_norm1(m_sub(A,B,B)) >= MACHEPS*A->m )
  413.         errmesg("mload()/m_save()");
  414.     }
  415.     }
  416.  
  417.     MEMCHK();
  418.  
  419.     /* Now, onto matrix factorisations */
  420.     A = m_resize(A,10,10);
  421.     B = m_resize(B,A->m,A->n);
  422.     m_copy(A,B);
  423.     x = v_resize(x,A->n);
  424.     y = v_resize(y,A->m);
  425.     z = v_resize(z,A->n);
  426.     u = v_resize(u,A->m);
  427.     v_rand(x);
  428.     mv_mlt(B,x,y);
  429.     z = v_copy(x,z);
  430.  
  431.     notice("LU factor/solve");
  432.     pivot = px_get(A->m);
  433.     LUfactor(A,pivot);
  434.     tracecatch(LUsolve(A,pivot,y,x),"main");
  435.     tracecatch(cond_est = LUcondest(A,pivot),"main");
  436.     printf("# cond(A) approx= %g\n", cond_est);
  437.     if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est)
  438.     {
  439.     errmesg("LUfactor()/LUsolve()");
  440.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  441.            v_norm2(v_sub(x,z,u)), MACHEPS);
  442.     }
  443.  
  444.     v_copy(y,x);
  445.     tracecatch(LUsolve(A,pivot,x,x),"main");
  446.     tracecatch(cond_est = LUcondest(A,pivot),"main");
  447.     if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est)
  448.     {
  449.     errmesg("LUfactor()/LUsolve()");
  450.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  451.            v_norm2(v_sub(x,z,u)), MACHEPS);
  452.     }
  453.  
  454.     vm_mlt(B,z,y);
  455.     v_copy(y,x);
  456.     tracecatch(LUTsolve(A,pivot,x,x),"main");
  457.     if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est)
  458.     {
  459.     errmesg("LUfactor()/LUTsolve()");
  460.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  461.            v_norm2(v_sub(x,z,u)), MACHEPS);
  462.     }
  463.  
  464.     MEMCHK();
  465.  
  466.     /* QR factorisation */
  467.     m_copy(B,A);
  468.     mv_mlt(B,z,y);
  469.     notice("QR factor/solve:");
  470.     diag = v_get(A->m);
  471.     beta = v_get(A->m);
  472.     QRfactor(A,diag);
  473.     QRsolve(A,diag,y,x);
  474.     if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est )
  475.     {
  476.     errmesg("QRfactor()/QRsolve()");
  477.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  478.            v_norm2(v_sub(x,z,u)), MACHEPS);
  479.     }
  480.     Q = m_get(A->m,A->m);
  481.     makeQ(A,diag,Q);
  482.     makeR(A,A);
  483.     m_mlt(Q,A,C);
  484.     m_sub(B,C,C);
  485.     if ( m_norm1(C) >= MACHEPS*m_norm1(Q)*m_norm1(B) )
  486.     {
  487.     errmesg("QRfactor()/makeQ()/makeR()");
  488.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  489.            m_norm1(C), MACHEPS);
  490.     }
  491.  
  492.     MEMCHK();
  493.  
  494.     /* now try with a non-square matrix */
  495.     A = m_resize(A,15,7);
  496.     m_rand(A);
  497.     B = m_copy(A,B);
  498.     diag = v_resize(diag,A->n);
  499.     beta = v_resize(beta,A->n);
  500.     x = v_resize(x,A->n);
  501.     y = v_resize(y,A->m);
  502.     v_rand(y);
  503.     QRfactor(A,diag);
  504.     x = QRsolve(A,diag,y,x);
  505.     /* z is the residual vector */
  506.     mv_mlt(B,x,z);    v_sub(z,y,z);
  507.     /* check B^T.z = 0 */
  508.     vm_mlt(B,z,u);
  509.     if ( v_norm2(u) >= MACHEPS*m_norm1(B)*v_norm2(y) )
  510.     {
  511.     errmesg("QRfactor()/QRsolve()");
  512.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  513.            v_norm2(u), MACHEPS);
  514.     }
  515.     Q = m_resize(Q,A->m,A->m);
  516.     makeQ(A,diag,Q);
  517.     makeR(A,A);
  518.     m_mlt(Q,A,C);
  519.     m_sub(B,C,C);
  520.     if ( m_norm1(C) >= MACHEPS*m_norm1(Q)*m_norm1(B) )
  521.     {
  522.     errmesg("QRfactor()/makeQ()/makeR()");
  523.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  524.            m_norm1(C), MACHEPS);
  525.     }
  526.     D = m_get(A->m,Q->m);
  527.     mtrm_mlt(Q,Q,D);
  528.     for ( i = 0; i < D->m; i++ )
  529.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  530.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q) )
  531.     {
  532.     errmesg("QRfactor()/makeQ()/makeR()");
  533.     printf("# QR orthogonality error = %g [cf MACHEPS = %g]\n",
  534.            m_norm1(D), MACHEPS);
  535.     }
  536.  
  537.     MEMCHK();
  538.  
  539.     /* QRCP factorisation */
  540.     m_copy(B,A);
  541.     notice("QR factor/solve with column pivoting");
  542.     pivot = px_resize(pivot,A->n);
  543.     QRCPfactor(A,diag,pivot);
  544.     z = v_resize(z,A->n);
  545.     QRCPsolve(A,diag,pivot,y,z);
  546.     /* pxinv_vec(pivot,z,x); */
  547.     /* now compute residual (z) vector */
  548.     mv_mlt(B,x,z);    v_sub(z,y,z);
  549.     /* check B^T.z = 0 */
  550.     vm_mlt(B,z,u);
  551.     if ( v_norm2(u) >= MACHEPS*m_norm1(B)*v_norm2(y) )
  552.     {
  553.     errmesg("QRCPfactor()/QRsolve()");
  554.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  555.            v_norm2(u), MACHEPS);
  556.     }
  557.  
  558.     Q = m_resize(Q,A->m,A->m);
  559.     makeQ(A,diag,Q);
  560.     makeR(A,A);
  561.     m_mlt(Q,A,C);
  562.     M_FREE(D);
  563.     D = m_get(B->m,B->n);
  564.     px_cols(pivot,C,D);
  565.     m_sub(B,D,D);
  566.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm1(B) )
  567.     {
  568.     errmesg("QRCPfactor()/makeQ()/makeR()");
  569.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  570.            m_norm1(D), MACHEPS);
  571.     }
  572.  
  573.     MEMCHK();
  574.  
  575.     /* Cholesky and LDL^T factorisation */
  576.     /* Use these for normal equations approach */
  577.     notice("Cholesky factor/solve");
  578.     mtrm_mlt(B,B,A);
  579.     CHfactor(A);
  580.     u = v_resize(u,B->n);
  581.     vm_mlt(B,y,u);
  582.     z = v_resize(z,B->n);
  583.     CHsolve(A,u,z);
  584.     v_sub(x,z,z);
  585.     if ( v_norm2(z) >= MACHEPS*v_norm2(x)*100 )
  586.     {
  587.     errmesg("CHfactor()/CHsolve()");
  588.     printf("# Cholesky solution error = %g [cf MACHEPS = %g]\n",
  589.            v_norm2(z), MACHEPS);
  590.     }
  591.     /* modified Cholesky factorisation should be identical with Cholesky
  592.        factorisation provided the matrix is "sufficiently positive definite */
  593.     mtrm_mlt(B,B,C);
  594.     MCHfactor(C,MACHEPS);
  595.     m_sub(A,C,C);
  596.     if ( m_norm1(C) >= MACHEPS*m_norm1(A) )
  597.     {
  598.     errmesg("MCHfactor()");
  599.     printf("# Modified Cholesky error = %g [cf MACHEPS = %g]\n",
  600.            m_norm1(C), MACHEPS);
  601.     }
  602.     /* now test the LDL^T factorisation -- using a negative def. matrix */
  603.     mtrm_mlt(B,B,A);
  604.     sm_mlt(-1.0,A,A);
  605.     m_copy(A,C);
  606.     LDLfactor(A);
  607.     LDLsolve(A,u,z);
  608.     w = v_get(A->m);
  609.     mv_mlt(C,z,w);
  610.     v_sub(w,u,w);
  611.     if ( v_norm2(w) >= MACHEPS*v_norm2(u)*m_norm1(C) )
  612.     {
  613.     errmesg("LDLfactor()/LDLsolve()");
  614.     printf("# LDL^T residual = %g [cf MACHEPS = %g]\n",
  615.            v_norm2(w), MACHEPS);
  616.     }
  617.     v_add(x,z,z);
  618.     if ( v_norm2(z) >= MACHEPS*v_norm2(x)*100 )
  619.     {
  620.     errmesg("LDLfactor()/LDLsolve()");
  621.     printf("# LDL^T solution error = %g [cf MACHEPS = %g]\n",
  622.            v_norm2(z), MACHEPS);
  623.     }
  624.  
  625.     MEMCHK();
  626.  
  627.     /* and now the Bunch-Kaufman-Parlett method */
  628.     /* set up D to be an indefinite diagonal matrix */
  629.     notice("Bunch-Kaufman-Parlett factor/solve");
  630.  
  631.     D = m_resize(D,B->m,B->m);
  632.     m_zero(D);
  633.     w = v_resize(w,B->m);
  634.     v_rand(w);
  635.     for ( i = 0; i < w->dim; i++ )
  636.     if ( v_entry(w,i) >= 0.5 )
  637.         m_set_val(D,i,i,1.0);
  638.     else
  639.         m_set_val(D,i,i,-1.0);
  640.     /* set A <- B^T.D.B */
  641.     C = m_resize(C,B->n,B->n);
  642.     C = mtrm_mlt(B,D,C);
  643.     A = m_mlt(C,B,A);
  644.     C = m_resize(C,B->n,B->n);
  645.     C = m_copy(A,C);
  646.     /* ... and use BKPfactor() */
  647.     blocks = px_get(A->m);
  648.     pivot = px_resize(pivot,A->m);
  649.     x = v_resize(x,A->m);
  650.     y = v_resize(y,A->m);
  651.     z = v_resize(z,A->m);
  652.     v_rand(x);
  653.     mv_mlt(A,x,y);
  654.     BKPfactor(A,pivot,blocks);
  655.     printf("# BKP pivot =\n");    px_output(pivot);
  656.     printf("# BKP blocks =\n");    px_output(blocks);
  657.     BKPsolve(A,pivot,blocks,y,z);
  658.     /* compute & check residual */
  659.     mv_mlt(C,z,w);
  660.     v_sub(w,y,w);
  661.     if ( v_norm2(w) >= MACHEPS*m_norm1(C)*v_norm2(z) )
  662.     {
  663.     errmesg("BKPfactor()/BKPsolve()");
  664.     printf("# BKP residual size = %g [cf MACHEPS = %g]\n",
  665.            v_norm2(w), MACHEPS);
  666.     }
  667.  
  668.     /* check update routines */
  669.     /* check LDLupdate() first */
  670.     notice("update L.D.L^T routine");
  671.     A = mtrm_mlt(B,B,A);
  672.     m_resize(C,A->m,A->n);
  673.     C = m_copy(A,C);
  674.     LDLfactor(A);
  675.     s1 = 3.7;
  676.     w = v_resize(w,A->m);
  677.     v_rand(w);
  678.     for ( i = 0; i < C->m; i++ )
  679.     for ( j = 0; j < C->n; j++ )
  680.         m_set_val(C,i,j,m_entry(C,i,j)+s1*v_entry(w,i)*v_entry(w,j));
  681.     LDLfactor(C);
  682.     LDLupdate(A,w,s1);
  683.     /* zero out strictly upper triangular parts of A and C */
  684.     for ( i = 0; i < A->m; i++ )
  685.     for ( j = i+1; j < A->n; j++ )
  686.     {
  687.         m_set_val(A,i,j,0.0);
  688.         m_set_val(C,i,j,0.0);
  689.     }
  690.     if ( m_norm1(m_sub(A,C,C)) >= sqrt(MACHEPS)*m_norm1(A) )
  691.     {
  692.     errmesg("LDLupdate()");
  693.     printf("# LDL update matrix error = %g [cf MACHEPS = %g]\n",
  694.            m_norm1(C), MACHEPS);
  695.     }
  696.  
  697.  
  698.     /* BAND MATRICES */
  699.  
  700. #define COL 40
  701. #define UDIAG  5
  702. #define LDIAG  2
  703.  
  704.    smrand(101);
  705.    bA = bd_get(LDIAG,UDIAG,COL);
  706.    bB = bd_get(LDIAG,UDIAG,COL);
  707.    bC = bd_get(LDIAG,UDIAG,COL);
  708.    A = m_resize(A,COL,COL);
  709.    B = m_resize(B,COL,COL);
  710.    pivot = px_resize(pivot,COL);
  711.    x = v_resize(x,COL);
  712.    w = v_resize(w,COL);
  713.    z = v_resize(z,COL);
  714.  
  715.    m_rand(A); 
  716.    /* generate band matrix */
  717.    mat2band(A,LDIAG,UDIAG,bA);
  718.    band2mat(bA,A);    /* now A is banded */
  719.    bB = bd_copy(bA,bB); 
  720.  
  721.    v_rand(x);  
  722.    mv_mlt(A,x,w);
  723.    z = v_copy(w,z);
  724.  
  725.    notice("band LU factorization");
  726.    bdLUfactor(bA,pivot);
  727.  
  728.    /* pivot will be changed */
  729.    bdLUsolve(bA,pivot,z,z);
  730.    v_sub(x,z,z);
  731.    if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
  732.       errmesg("incorrect solution (band LU factorization)");
  733.       printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
  734.          v_norm2(z),MACHEPS);
  735.    }
  736.  
  737.    /* solve transpose system */
  738.  
  739.    notice("band LU factorization for transpose system");
  740.    m_transp(A,B);
  741.    mv_mlt(B,x,w);
  742.  
  743.    bd_copy(bB,bA);
  744.    bd_transp(bA,bA);  
  745.    /* transposition in situ */
  746.    bd_transp(bA,bB);
  747.    bd_transp(bB,bB);
  748.  
  749.    bdLUfactor(bB,pivot);
  750.  
  751.    bdLUsolve(bB,pivot,w,z);
  752.    v_sub(x,z,z);
  753.    if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
  754.       errmesg("incorrect solution (band transposed LU factorization)");
  755.       printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
  756.          v_norm2(z),MACHEPS);
  757.    }
  758.  
  759.  
  760.    /* Cholesky factorization */
  761.  
  762.    notice("band Choleski LDL' factorization");
  763.    m_add(A,B,A);  /* symmetric matrix */
  764.    for (i=0; i < COL; i++)     /* positive definite */
  765.      A->me[i][i] += 2*LDIAG;   
  766.  
  767.    mat2band(A,LDIAG,LDIAG,bA);
  768.    band2mat(bA,A);              /* corresponding matrix A */
  769.  
  770.    v_rand(x);
  771.    mv_mlt(A,x,w);
  772.    z = v_copy(w,z);
  773.    
  774.    bdLDLfactor(bA);
  775.  
  776.    z = bdLDLsolve(bA,z,z);
  777.    v_sub(x,z,z);
  778.    if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
  779.       errmesg("incorrect solution (band LDL' factorization)");
  780.       printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
  781.          v_norm2(z),MACHEPS);
  782.    }
  783.  
  784.    /* new bandwidths */
  785.    m_rand(A);
  786.    bA = bd_resize(bA,UDIAG,LDIAG,COL);
  787.    bB = bd_resize(bB,UDIAG,LDIAG,COL);
  788.    mat2band(A,UDIAG,LDIAG,bA);
  789.    band2mat(bA,A);
  790.    bd_copy(bA,bB);
  791.  
  792.    mv_mlt(A,x,w);
  793.  
  794.    notice("band LU factorization (resized)");
  795.    bdLUfactor(bA,pivot);
  796.  
  797.    /* pivot will be changed */
  798.    bdLUsolve(bA,pivot,w,z);
  799.    v_sub(x,z,z);
  800.    if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
  801.       errmesg("incorrect solution (band LU factorization)");
  802.       printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
  803.          v_norm2(z),MACHEPS);
  804.    }
  805.  
  806.    /* testing transposition */
  807.  
  808.    notice("band matrix transposition");
  809.    m_zero(bA->mat);
  810.    bd_copy(bB,bA);
  811.    m_zero(bB->mat);
  812.    bd_copy(bA,bB);
  813.  
  814.    bd_transp(bB,bB);
  815.    bd_transp(bB,bB);
  816.  
  817.    m_zero(bC->mat);
  818.    bd_copy(bB,bC);
  819.  
  820.    m_sub(bA->mat,bC->mat,bC->mat);
  821.    if (m_norm_inf(bC->mat) > MACHEPS*bC->mat->n) {
  822.       errmesg("band transposition");
  823.       printf(" difference ||A - (A')'|| = %g\n",m_norm_inf(bC->mat));
  824.    }
  825.  
  826.    bd_free(bA);
  827.    bd_free(bB);
  828.    bd_free(bC);
  829.  
  830.  
  831.     MEMCHK();
  832.  
  833.     /* now check QRupdate() routine */
  834.     notice("update QR routine");
  835.  
  836.     B = m_resize(B,15,7);
  837.     A = m_resize(A,B->m,B->n);
  838.     m_copy(B,A);
  839.     diag = v_resize(diag,A->n);
  840.     beta = v_resize(beta,A->n);
  841.     QRfactor(A,diag);
  842.     Q = m_resize(Q,A->m,A->m);
  843.     makeQ(A,diag,Q);
  844.     makeR(A,A);
  845.     m_resize(C,A->m,A->n);
  846.     w = v_resize(w,A->m);
  847.     v = v_resize(v,A->n);
  848.     u = v_resize(u,A->m);
  849.     v_rand(w);
  850.     v_rand(v);
  851.     vm_mlt(Q,w,u);
  852.     QRupdate(Q,A,u,v);
  853.     m_mlt(Q,A,C);
  854.     for ( i = 0; i < B->m; i++ )
  855.     for ( j = 0; j < B->n; j++ )
  856.         m_set_val(B,i,j,m_entry(B,i,j)+v_entry(w,i)*v_entry(v,j));
  857.     m_sub(B,C,C);
  858.     if ( m_norm1(C) >= MACHEPS*m_norm1(A)*m_norm1(Q)*2 )
  859.     {
  860.     errmesg("QRupdate()");
  861.     printf("# Reconstruction error in QR update = %g [cf MACHEPS = %g]\n",
  862.            m_norm1(C), MACHEPS);
  863.     }
  864.     m_resize(D,Q->m,Q->n);
  865.     mtrm_mlt(Q,Q,D);
  866.     for ( i = 0; i < D->m; i++ )
  867.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  868.     if ( m_norm1(D) >= 10*MACHEPS*m_norm1(Q)*m_norm_inf(Q) )
  869.     {
  870.     errmesg("QRupdate()");
  871.     printf("# QR update orthogonality error = %g [cf MACHEPS = %g]\n",
  872.            m_norm1(D), MACHEPS);
  873.     }
  874.  
  875.     /* Now check eigenvalue/SVD routines */
  876.     notice("eigenvalue and SVD routines");
  877.     A = m_resize(A,11,11);
  878.     B = m_resize(B,A->m,A->n);
  879.     C = m_resize(C,A->m,A->n);
  880.     D = m_resize(D,A->m,A->n);
  881.     Q = m_resize(Q,A->m,A->n);
  882.  
  883.     m_rand(A);
  884.     /* A <- A + A^T  for symmetric case */
  885.     m_add(A,m_transp(A,C),A);
  886.     u = v_resize(u,A->m);
  887.     u = symmeig(A,Q,u);
  888.     m_zero(B);
  889.     for ( i = 0; i < B->m; i++ )
  890.     m_set_val(B,i,i,v_entry(u,i));
  891.     m_mlt(Q,B,C);
  892.     mmtr_mlt(C,Q,D);
  893.     m_sub(A,D,D);
  894.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*v_norm_inf(u)*3 )
  895.     {
  896.     errmesg("symmeig()");
  897.     printf("# Reconstruction error = %g [cf MACHEPS = %g]\n",
  898.            m_norm1(D), MACHEPS);
  899.     }
  900.     mtrm_mlt(Q,Q,D);
  901.     for ( i = 0; i < D->m; i++ )
  902.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  903.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*3 )
  904.     {
  905.     errmesg("symmeig()");
  906.     printf("# symmeig() orthogonality error = %g [cf MACHEPS = %g]\n",
  907.            m_norm1(D), MACHEPS);
  908.     }
  909.  
  910.     MEMCHK();
  911.  
  912.     /* now test (real) Schur decomposition */
  913.     /* m_copy(A,B); */
  914.     M_FREE(A);
  915.     A = m_get(11,11);
  916.     m_rand(A);
  917.     B = m_copy(A,B);
  918.     MEMCHK();
  919.  
  920.     B = schur(B,Q);
  921.     MEMCHK();
  922.  
  923.     m_mlt(Q,B,C);
  924.     mmtr_mlt(C,Q,D);
  925.     MEMCHK();
  926.     m_sub(A,D,D);
  927.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*m_norm1(B)*5 )
  928.     {
  929.     errmesg("schur()");
  930.     printf("# Schur reconstruction error = %g [cf MACHEPS = %g]\n",
  931.            m_norm1(D), MACHEPS);
  932.     }
  933.  
  934.     /* orthogonality check */
  935.     mmtr_mlt(Q,Q,D);
  936.     for ( i = 0; i < D->m; i++ )
  937.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  938.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*10 )
  939.     {
  940.     errmesg("schur()");
  941.     printf("# Schur orthogonality error = %g [cf MACHEPS = %g]\n",
  942.            m_norm1(D), MACHEPS);
  943.     }
  944.  
  945.     MEMCHK();
  946.  
  947.     /* now test SVD */
  948.     A = m_resize(A,11,7);
  949.     m_rand(A);
  950.     U = m_get(A->n,A->n);
  951.     Q = m_resize(Q,A->m,A->m);
  952.     u = v_resize(u,max(A->m,A->n));
  953.     svd(A,Q,U,u);
  954.     /* check reconstruction of A */
  955.     D = m_resize(D,A->m,A->n);
  956.     C = m_resize(C,A->m,A->n);
  957.     m_zero(D);
  958.     for ( i = 0; i < min(A->m,A->n); i++ )
  959.     m_set_val(D,i,i,v_entry(u,i));
  960.     mtrm_mlt(Q,D,C);
  961.     m_mlt(C,U,D);
  962.     m_sub(A,D,D);
  963.     if ( m_norm1(D) >= MACHEPS*m_norm1(U)*m_norm_inf(Q)*m_norm1(A) )
  964.     {
  965.     errmesg("svd()");
  966.     printf("# SVD reconstruction error = %g [cf MACHEPS = %g]\n",
  967.            m_norm1(D), MACHEPS);
  968.     }
  969.     /* check orthogonality of Q and U */
  970.     D = m_resize(D,Q->n,Q->n);
  971.     mtrm_mlt(Q,Q,D);
  972.     for ( i = 0; i < D->m; i++ )
  973.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  974.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*5 )
  975.     {
  976.     errmesg("svd()");
  977.     printf("# SVD orthognality error (Q) = %g [cf MACHEPS = %g\n",
  978.            m_norm1(D), MACHEPS);
  979.     }
  980.     D = m_resize(D,U->n,U->n);
  981.     mtrm_mlt(U,U,D);
  982.     for ( i = 0; i < D->m; i++ )
  983.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  984.     if ( m_norm1(D) >= MACHEPS*m_norm1(U)*m_norm_inf(U)*5 )
  985.     {
  986.     errmesg("svd()");
  987.     printf("# SVD orthognality error (U) = %g [cf MACHEPS = %g\n",
  988.            m_norm1(D), MACHEPS);
  989.     }
  990.     for ( i = 0; i < u->dim; i++ )
  991.     if ( v_entry(u,i) < 0 || (i < u->dim-1 &&
  992.                   v_entry(u,i+1) > v_entry(u,i)) )
  993.         break;
  994.     if ( i < u->dim )
  995.     {
  996.     errmesg("svd()");
  997.     printf("# SVD sorting error\n");
  998.     }
  999.  
  1000.  
  1001.     /* test of long vectors */
  1002.     notice("Long vectors");
  1003.     x = v_resize(x,100000);
  1004.     y = v_resize(y,100000);
  1005.     z = v_resize(z,100000);
  1006.     v_rand(x);
  1007.     v_rand(y);
  1008.     v_mltadd(x,y,3.0,z);
  1009.     sv_mlt(1.0/3.0,z,z);
  1010.     v_mltadd(z,x,-1.0/3.0,z);
  1011.     v_sub(z,y,x);
  1012.     if (v_norm2(x) >= MACHEPS*(x->dim)) {
  1013.        errmesg("long vectors");
  1014.        printf(" norm = %g\n",v_norm2(x));
  1015.     }
  1016.  
  1017.     mem_stat_free(1);
  1018.  
  1019.     MEMCHK();
  1020.  
  1021.     /**************************************************
  1022.     VEC        *x, *y, *z, *u, *v, *w;
  1023.     VEC        *diag, *beta;
  1024.     PERM    *pi1, *pi2, *pi3, *pivot, *blocks;
  1025.     MAT        *A, *B, *C, *D, *Q, *U;
  1026.     **************************************************/
  1027.     V_FREE(x);        V_FREE(y);    V_FREE(z);
  1028.     V_FREE(u);        V_FREE(v);    V_FREE(w);
  1029.     V_FREE(diag);    V_FREE(beta);
  1030.     PX_FREE(pi1);    PX_FREE(pi2);    PX_FREE(pi3);
  1031.     PX_FREE(pivot);    PX_FREE(blocks);
  1032.     M_FREE(A);        M_FREE(B);    M_FREE(C);
  1033.     M_FREE(D);        M_FREE(Q);    M_FREE(U);
  1034.  
  1035.     MEMCHK();
  1036.     printf("# Finished torture test\n");
  1037.     mem_info();
  1038.  
  1039.     return 0;
  1040. }
  1041.  
  1042.