home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / sptort.c < prev    next >
C/C++ Source or Header  |  1994-01-14  |  12KB  |  484 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 tests for the sparse matrix part of Meschach
  29. */
  30.  
  31. #include    <stdio.h>
  32. #include    <math.h>
  33. #include    "matrix2.h"
  34. #include    "sparse2.h"
  35. #include        "iter.h"
  36.  
  37. #define    errmesg(mesg)    printf("Error: %s error: line %d\n",mesg,__LINE__)
  38. #define notice(mesg)    printf("# Testing %s...\n",mesg);
  39.  
  40. /* for iterative methods */
  41.  
  42. #if REAL == DOUBLE
  43. #define    EPS    1e-7
  44. #elif REAL == FLOAT
  45. #define EPS   1e-3
  46. #endif
  47.  
  48. int    chk_col_access(A)
  49. SPMAT    *A;
  50. {
  51.     int        i, j, nxt_idx, nxt_row, scan_cnt, total_cnt;
  52.     SPROW    *r;
  53.     row_elt    *e;
  54.  
  55.     if ( ! A )
  56.     error(E_NULL,"chk_col_access");
  57.     if ( ! A->flag_col )
  58.     return FALSE;
  59.  
  60.     /* scan down each column, counting the number of entries met */
  61.     scan_cnt = 0;
  62.     for ( j = 0; j < A->n; j++ )
  63.     {
  64.     i = -1;
  65.     nxt_idx = A->start_idx[j];
  66.     nxt_row = A->start_row[j];
  67.     while ( nxt_row >= 0 && nxt_idx >= 0 && nxt_row > i )
  68.     {
  69.         i = nxt_row;
  70.         r = &(A->row[i]);
  71.         e = &(r->elt[nxt_idx]);
  72.         nxt_idx = e->nxt_idx;
  73.         nxt_row = e->nxt_row;
  74.         scan_cnt++;
  75.     }
  76.     }
  77.  
  78.     total_cnt = 0;
  79.     for ( i = 0; i < A->m; i++ )
  80.     total_cnt += A->row[i].len;
  81.     if ( total_cnt != scan_cnt )
  82.     return FALSE;
  83.     else
  84.     return TRUE;
  85. }
  86.  
  87.  
  88. void    main(argc, argv)
  89. int    argc;
  90. char    *argv[];
  91. {
  92.     VEC        *x, *y, *z, *u, *v;
  93.     Real    s1, s2;
  94.     PERM    *pivot;
  95.     SPMAT    *A, *B, *C;
  96.     SPMAT       *B1, *C1;
  97.     SPROW    *r;
  98.     int        i, j, k, deg, seed, m, n;
  99.  
  100.  
  101.     mem_info_on(TRUE);
  102.  
  103.     setbuf(stdout, (char *)NULL);
  104.     /* get seed if in argument list */
  105.     if ( argc == 1 )
  106.     seed = 1111;
  107.     else if ( argc == 2 && sscanf(argv[1],"%d",&seed) == 1 )
  108.     ;
  109.     else
  110.     {
  111.     printf("usage: %s [seed]\n", argv[0]);
  112.     exit(0);
  113.     }
  114.     srand(seed);
  115.  
  116.     /* set up two random sparse matrices */
  117.     m = 120;
  118.     n = 100;
  119.     deg = 8;
  120.     notice("allocating sparse matrices");
  121.     A = sp_get(m,n,deg);
  122.     B = sp_get(m,n,deg);
  123.     notice("setting and getting matrix entries");
  124.     for ( k = 0; k < m*deg; k++ )
  125.     {
  126.     i = (rand() >> 8) % m;
  127.     j = (rand() >> 8) % n;
  128.     sp_set_val(A,i,j,rand()/((Real)MAX_RAND));
  129.     i = (rand() >> 8) % m;
  130.     j = (rand() >> 8) % n;
  131.     sp_set_val(B,i,j,rand()/((Real)MAX_RAND));
  132.     }
  133.     for ( k = 0; k < 10; k++ )
  134.     {
  135.     s1 = rand()/((Real)MAX_RAND);
  136.     i = (rand() >> 8) % m;
  137.     j = (rand() >> 8) % n;
  138.     sp_set_val(A,i,j,s1);
  139.     s2 = sp_get_val(A,i,j);
  140.     if ( fabs(s1 - s2) >= MACHEPS )
  141.         break;
  142.     }
  143.     if ( k < 10 )
  144.     errmesg("sp_set_val()/sp_get_val()");
  145.  
  146.     /* test copy routines */
  147.     notice("copy routines");
  148.     x = v_get(n);
  149.     y = v_get(m);
  150.     z = v_get(m);
  151.     /* first copy routine */
  152.     C = sp_copy(A);
  153.     for ( k = 0; k < 100; k++ )
  154.     {
  155.     v_rand(x);
  156.     sp_mv_mlt(A,x,y);
  157.     sp_mv_mlt(C,x,z);
  158.     if ( v_norm_inf(v_sub(y,z,z)) >= MACHEPS*deg*m )
  159.         break;
  160.     }
  161.     if ( k < 100 )
  162.     {
  163.     errmesg("sp_copy()/sp_mv_mlt()");
  164.     printf("# Error in A.x (inf norm) = %g [cf MACHEPS = %g]\n",
  165.            v_norm_inf(z), MACHEPS);
  166.     }
  167.     /* second copy routine
  168.        -- note that A & B have different sparsity patterns */
  169.  
  170.     mem_stat_mark(1);
  171.     sp_copy2(A,B);
  172.     mem_stat_free(1);
  173.     for ( k = 0; k < 10; k++ )
  174.     {
  175.     v_rand(x);
  176.     sp_mv_mlt(A,x,y);
  177.     sp_mv_mlt(B,x,z);
  178.     if ( v_norm_inf(v_sub(y,z,z)) >= MACHEPS*deg*m )
  179.         break;
  180.     }
  181.     if ( k < 10 )
  182.     {
  183.     errmesg("sp_copy2()/sp_mv_mlt()");
  184.     printf("# Error in A.x (inf norm) = %g [cf MACHEPS = %g]\n",
  185.            v_norm_inf(z), MACHEPS);
  186.     }
  187.  
  188.     /* now check compacting routine */
  189.     notice("compacting routine");
  190.     sp_compact(B,0.0);
  191.     for ( k = 0; k < 10; k++ )
  192.     {
  193.     v_rand(x);
  194.     sp_mv_mlt(A,x,y);
  195.     sp_mv_mlt(B,x,z);
  196.     if ( v_norm_inf(v_sub(y,z,z)) >= MACHEPS*deg*m )
  197.         break;
  198.     }
  199.     if ( k < 10 )
  200.     {
  201.     errmesg("sp_compact()");
  202.     printf("# Error in A.x (inf norm) = %g [cf MACHEPS = %g]\n",
  203.            v_norm_inf(z), MACHEPS);
  204.     }
  205.     for ( i = 0; i < B->m; i++ )
  206.     {
  207.     r = &(B->row[i]);
  208.     for ( j = 0; j < r->len; j++ )
  209.         if ( r->elt[j].val == 0.0 )
  210.         break;
  211.     }
  212.     if ( i < B->m )
  213.     {
  214.     errmesg("sp_compact()");
  215.     printf("# Zero entry in compacted matrix\n");
  216.     }
  217.  
  218.     /* check column access paths */
  219.     notice("resizing and access paths");
  220.     sp_resize(A,A->m+10,A->n+10);
  221.     for ( k = 0 ; k < 20; k++ )
  222.     {
  223.     i = A->m - 1 - ((rand() >> 8) % 10);
  224.     j = A->n - 1 - ((rand() >> 8) % 10);
  225.     s1 = rand()/((Real)MAX_RAND);
  226.     sp_set_val(A,i,j,s1);
  227.     if ( fabs(s1 - sp_get_val(A,i,j)) >= MACHEPS )
  228.         break;
  229.     }
  230.     if ( k < 20 )
  231.     errmesg("sp_resize()");
  232.     sp_col_access(A);
  233.     if ( ! chk_col_access(A) )
  234.     {
  235.     errmesg("sp_col_access()");
  236.     }
  237.     sp_diag_access(A);
  238.     for ( i = 0; i < A->m; i++ )
  239.     {
  240.     r = &(A->row[i]);
  241.     if ( r->diag != sprow_idx(r,i) )
  242.         break;
  243.     }
  244.     if ( i < A->m )
  245.     {
  246.     errmesg("sp_diag_access()");
  247.     }
  248.  
  249.     /* test both sp_mv_mlt() and sp_vm_mlt() */
  250.     x = v_resize(x,B->n);
  251.     y = v_resize(y,B->m);
  252.     u = v_get(B->m);
  253.     v = v_get(B->n);
  254.     for ( k = 0; k < 10; k++ )
  255.     {
  256.     v_rand(x);
  257.     v_rand(y);
  258.     sp_mv_mlt(B,x,u);
  259.     sp_vm_mlt(B,y,v);
  260.     if ( fabs(in_prod(x,v) - in_prod(y,u)) >=
  261.         MACHEPS*v_norm2(x)*v_norm2(u)*5 )
  262.         break;
  263.     }
  264.     if ( k < 10 )
  265.     {
  266.     errmesg("sp_mv_mlt()/sp_vm_mlt()");
  267.     printf("# Error in inner products = %g [cf MACHEPS = %g]\n",
  268.            fabs(in_prod(x,v) - in_prod(y,u)), MACHEPS);
  269.     }
  270.  
  271.     SP_FREE(A);
  272.     SP_FREE(B);
  273.     SP_FREE(C);
  274.  
  275.     /* now test Cholesky and LU factorise and solve */
  276.     notice("sparse Cholesky factorise/solve");
  277.     A = iter_gen_sym(120,8);
  278.     B = sp_copy(A);
  279.     spCHfactor(A);
  280.     x = v_resize(x,A->m);
  281.     y = v_resize(y,A->m);
  282.     v_rand(x);
  283.     sp_mv_mlt(B,x,y);
  284.     z = v_resize(z,A->m);
  285.     spCHsolve(A,y,z);
  286.     v = v_resize(v,A->m);
  287.     sp_mv_mlt(B,z,v);
  288.     /* compute residual */
  289.     v_sub(y,v,v);
  290.     if ( v_norm2(v) >= MACHEPS*v_norm2(y)*10 )
  291.     {
  292.     errmesg("spCHfactor()/spCHsolve()");
  293.     printf("# Sparse Cholesky residual = %g [cf MACHEPS = %g]\n",
  294.            v_norm2(v), MACHEPS);
  295.     }
  296.     /* compute error in solution */
  297.     v_sub(x,z,z);
  298.     if ( v_norm2(z) > MACHEPS*v_norm2(x)*10 )
  299.     {
  300.     errmesg("spCHfactor()/spCHsolve()");
  301.     printf("# Solution error = %g [cf MACHEPS = %g]\n",
  302.            v_norm2(z), MACHEPS);
  303.     }
  304.  
  305.     /* now test symbolic and incomplete factorisation */
  306.     SP_FREE(A);
  307.     A = sp_copy(B);
  308.     
  309.     mem_stat_mark(2);
  310.     spCHsymb(A);
  311.     mem_stat_mark(2);
  312.  
  313.     spICHfactor(A);
  314.     spCHsolve(A,y,z);
  315.     v = v_resize(v,A->m);
  316.     sp_mv_mlt(B,z,v);
  317.     /* compute residual */
  318.     v_sub(y,v,v);
  319.     if ( v_norm2(v) >= MACHEPS*v_norm2(y)*5 )
  320.     {
  321.     errmesg("spCHsymb()/spICHfactor()");
  322.     printf("# Sparse Cholesky residual = %g [cf MACHEPS = %g]\n",
  323.            v_norm2(v), MACHEPS);
  324.     }
  325.     /* compute error in solution */
  326.     v_sub(x,z,z);
  327.     if ( v_norm2(z) > MACHEPS*v_norm2(x)*10 )
  328.     {
  329.     errmesg("spCHsymb()/spICHfactor()");
  330.     printf("# Solution error = %g [cf MACHEPS = %g]\n",
  331.            v_norm2(z), MACHEPS);
  332.     }
  333.  
  334.     /* now test sparse LU factorisation */
  335.     notice("sparse LU factorise/solve");
  336.     SP_FREE(A);
  337.     SP_FREE(B);
  338.     A = iter_gen_nonsym(100,100,8,1.0);
  339.  
  340.     B = sp_copy(A);
  341.     x = v_resize(x,A->n);
  342.     z = v_resize(z,A->n);
  343.     y = v_resize(y,A->m);
  344.     v = v_resize(v,A->m);
  345.  
  346.     v_rand(x);
  347.     sp_mv_mlt(B,x,y);
  348.     pivot = px_get(A->m);
  349.  
  350.     mem_stat_mark(3);
  351.     spLUfactor(A,pivot,0.1);
  352.     spLUsolve(A,pivot,y,z);
  353.     mem_stat_free(3);
  354.     sp_mv_mlt(B,z,v);
  355.  
  356.     /* compute residual */
  357.     v_sub(y,v,v);
  358.     if ( v_norm2(v) >= MACHEPS*v_norm2(y)*A->m )
  359.     {
  360.     errmesg("spLUfactor()/spLUsolve()");
  361.     printf("# Sparse LU residual = %g [cf MACHEPS = %g]\n",
  362.            v_norm2(v), MACHEPS);
  363.     }
  364.     /* compute error in solution */
  365.     v_sub(x,z,z);
  366.     if ( v_norm2(z) > MACHEPS*v_norm2(x)*100*A->m )
  367.     {
  368.     errmesg("spLUfactor()/spLUsolve()");
  369.     printf("# Sparse LU solution error = %g [cf MACHEPS = %g]\n",
  370.            v_norm2(z), MACHEPS);
  371.     }
  372.  
  373.     /* now check spLUTsolve */
  374.     mem_stat_mark(4);
  375.     sp_vm_mlt(B,x,y);
  376.     spLUTsolve(A,pivot,y,z);
  377.     sp_vm_mlt(B,z,v);
  378.     mem_stat_free(4);
  379.  
  380.     /* compute residual */
  381.     v_sub(y,v,v);
  382.     if ( v_norm2(v) >= MACHEPS*v_norm2(y)*A->m )
  383.     {
  384.     errmesg("spLUTsolve()");
  385.     printf("# Sparse LU residual = %g [cf MACHEPS = %g]\n",
  386.            v_norm2(v), MACHEPS);
  387.     }
  388.     /* compute error in solution */
  389.     v_sub(x,z,z);
  390.     if ( v_norm2(z) > MACHEPS*v_norm2(x)*100*A->m )
  391.     {
  392.     errmesg("spLUTsolve()");
  393.     printf("# Sparse LU solution error = %g [cf MACHEPS = %g]\n",
  394.            v_norm2(z), MACHEPS);
  395.     }
  396.  
  397.     /* algebraic operations */
  398.     notice("addition,subtraction and multiplying by a number");
  399.     SP_FREE(A);
  400.     SP_FREE(B);
  401.  
  402.     m = 120;
  403.     n = 120;
  404.     deg = 5;
  405.     A = sp_get(m,n,deg);
  406.     B = sp_get(m,n,deg);
  407.     C = sp_get(m,n,deg);
  408.     C1 = sp_get(m,n,deg);
  409.  
  410.     for ( k = 0; k < m*deg; k++ )
  411.     {
  412.     i = (rand() >> 8) % m;
  413.     j = (rand() >> 8) % n;
  414.     sp_set_val(A,i,j,rand()/((Real)MAX_RAND));
  415.     i = (rand() >> 8) % m;
  416.     j = (rand() >> 8) % n;
  417.     sp_set_val(B,i,j,rand()/((Real)MAX_RAND));
  418.     }
  419.     
  420.     s1 = mrand(); 
  421.     B1 = sp_copy(B);
  422.  
  423.     mem_stat_mark(1);
  424.     sp_smlt(B,s1,C);
  425.     sp_add(A,C,C1);
  426.     sp_sub(C1,A,C);
  427.     sp_smlt(C,-1.0/s1,C);
  428.     sp_add(C,B1,C);
  429.  
  430.     s2 = 0.0;
  431.     for (k=0; k < C->m; k++) {
  432.        r = &(C->row[i]);
  433.        for (j=0; j < r->len; j++) {
  434.       if (s2 < fabs(r->elt[j].val)) 
  435.         s2 = fabs(r->elt[j].val);
  436.        }
  437.     }
  438.  
  439.     if (s2 > MACHEPS*A->m) {
  440.        errmesg("add, sub, mlt sparse matrices (args not in situ)\n");
  441.        printf(" difference = %g [MACEPS = %g]\n",s2,MACHEPS);
  442.     }
  443.  
  444.     sp_mltadd(A,B1,s1,C1);
  445.     sp_sub(C1,A,A);
  446.     sp_smlt(A,1.0/s1,C1);
  447.     sp_sub(C1,B1,C1);
  448.     mem_stat_free(1);
  449.  
  450.     s2 = 0.0;
  451.     for (k=0; k < C1->m; k++) {
  452.        r = &(C1->row[i]);
  453.        for (j=0; j < r->len; j++) {
  454.       if (s2 < fabs(r->elt[j].val)) 
  455.         s2 = fabs(r->elt[j].val);
  456.        }
  457.     }
  458.  
  459.     if (s2 > MACHEPS*A->m) {
  460.        errmesg("add, sub, mlt sparse matrices (args not in situ)\n");
  461.        printf(" difference = %g [MACEPS = %g]\n",s2,MACHEPS);
  462.     }
  463.  
  464.     V_FREE(x);
  465.     V_FREE(y);    
  466.     V_FREE(z);
  467.     V_FREE(u);
  468.     V_FREE(v);  
  469.     PX_FREE(pivot);
  470.     SP_FREE(A);
  471.     SP_FREE(B);
  472.     SP_FREE(C);
  473.     SP_FREE(B1);
  474.     SP_FREE(C1);
  475.  
  476.     printf("# Done testing (%s)\n",argv[0]);
  477.     mem_info();
  478. }
  479.     
  480.  
  481.  
  482.  
  483.  
  484.