home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / sptort < prev    next >
Encoding:
Text File  |  1994-02-28  |  11.0 KB  |  486 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, m_old, n, n_old;
  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.     m_old = A->m-1;
  221.     n_old = A->n-1;
  222.     A = sp_resize(A,A->m+10,A->n+10);
  223.     for ( k = 0 ; k < 20; k++ )
  224.     {
  225.     i = m_old + ((rand() >> 8) % 10);
  226.     j = n_old + ((rand() >> 8) % 10);
  227.     s1 = rand()/((Real)MAX_RAND);
  228.     sp_set_val(A,i,j,s1);
  229.     if ( fabs(s1 - sp_get_val(A,i,j)) >= MACHEPS )
  230.         break;
  231.     }
  232.     if ( k < 20 )
  233.     errmesg("sp_resize()");
  234.     sp_col_access(A);
  235.     if ( ! chk_col_access(A) )
  236.     {
  237.     errmesg("sp_col_access()");
  238.     }
  239.     sp_diag_access(A);
  240.     for ( i = 0; i < A->m; i++ )
  241.     {
  242.     r = &(A->row[i]);
  243.     if ( r->diag != sprow_idx(r,i) )
  244.         break;
  245.     }
  246.     if ( i < A->m )
  247.     {
  248.     errmesg("sp_diag_access()");
  249.     }
  250.  
  251.     /* test both sp_mv_mlt() and sp_vm_mlt() */
  252.     x = v_resize(x,B->n);
  253.     y = v_resize(y,B->m);
  254.     u = v_get(B->m);
  255.     v = v_get(B->n);
  256.     for ( k = 0; k < 10; k++ )
  257.     {
  258.     v_rand(x);
  259.     v_rand(y);
  260.     sp_mv_mlt(B,x,u);
  261.     sp_vm_mlt(B,y,v);
  262.     if ( fabs(in_prod(x,v) - in_prod(y,u)) >=
  263.         MACHEPS*v_norm2(x)*v_norm2(u)*5 )
  264.         break;
  265.     }
  266.     if ( k < 10 )
  267.     {
  268.     errmesg("sp_mv_mlt()/sp_vm_mlt()");
  269.     printf("# Error in inner products = %g [cf MACHEPS = %g]\n",
  270.            fabs(in_prod(x,v) - in_prod(y,u)), MACHEPS);
  271.     }
  272.  
  273.     SP_FREE(A);
  274.     SP_FREE(B);
  275.     SP_FREE(C);
  276.  
  277.     /* now test Cholesky and LU factorise and solve */
  278.     notice("sparse Cholesky factorise/solve");
  279.     A = iter_gen_sym(120,8);
  280.     B = sp_copy(A);
  281.     spCHfactor(A);
  282.     x = v_resize(x,A->m);
  283.     y = v_resize(y,A->m);
  284.     v_rand(x);
  285.     sp_mv_mlt(B,x,y);
  286.     z = v_resize(z,A->m);
  287.     spCHsolve(A,y,z);
  288.     v = v_resize(v,A->m);
  289.     sp_mv_mlt(B,z,v);
  290.     /* compute residual */
  291.     v_sub(y,v,v);
  292.     if ( v_norm2(v) >= MACHEPS*v_norm2(y)*10 )
  293.     {
  294.     errmesg("spCHfactor()/spCHsolve()");
  295.     printf("# Sparse Cholesky residual = %g [cf MACHEPS = %g]\n",
  296.            v_norm2(v), MACHEPS);
  297.     }
  298.     /* compute error in solution */
  299.     v_sub(x,z,z);
  300.     if ( v_norm2(z) > MACHEPS*v_norm2(x)*10 )
  301.     {
  302.     errmesg("spCHfactor()/spCHsolve()");
  303.     printf("# Solution error = %g [cf MACHEPS = %g]\n",
  304.            v_norm2(z), MACHEPS);
  305.     }
  306.  
  307.     /* now test symbolic and incomplete factorisation */
  308.     SP_FREE(A);
  309.     A = sp_copy(B);
  310.     
  311.     mem_stat_mark(2);
  312.     spCHsymb(A);
  313.     mem_stat_mark(2);
  314.  
  315.     spICHfactor(A);
  316.     spCHsolve(A,y,z);
  317.     v = v_resize(v,A->m);
  318.     sp_mv_mlt(B,z,v);
  319.     /* compute residual */
  320.     v_sub(y,v,v);
  321.     if ( v_norm2(v) >= MACHEPS*v_norm2(y)*5 )
  322.     {
  323.     errmesg("spCHsymb()/spICHfactor()");
  324.     printf("# Sparse Cholesky residual = %g [cf MACHEPS = %g]\n",
  325.            v_norm2(v), MACHEPS);
  326.     }
  327.     /* compute error in solution */
  328.     v_sub(x,z,z);
  329.     if ( v_norm2(z) > MACHEPS*v_norm2(x)*10 )
  330.     {
  331.     errmesg("spCHsymb()/spICHfactor()");
  332.     printf("# Solution error = %g [cf MACHEPS = %g]\n",
  333.            v_norm2(z), MACHEPS);
  334.     }
  335.  
  336.     /* now test sparse LU factorisation */
  337.     notice("sparse LU factorise/solve");
  338.     SP_FREE(A);
  339.     SP_FREE(B);
  340.     A = iter_gen_nonsym(100,100,8,1.0);
  341.  
  342.     B = sp_copy(A);
  343.     x = v_resize(x,A->n);
  344.     z = v_resize(z,A->n);
  345.     y = v_resize(y,A->m);
  346.     v = v_resize(v,A->m);
  347.  
  348.     v_rand(x);
  349.     sp_mv_mlt(B,x,y);
  350.     pivot = px_get(A->m);
  351.  
  352.     mem_stat_mark(3);
  353.     spLUfactor(A,pivot,0.1);
  354.     spLUsolve(A,pivot,y,z);
  355.     mem_stat_free(3);
  356.     sp_mv_mlt(B,z,v);
  357.  
  358.     /* compute residual */
  359.     v_sub(y,v,v);
  360.     if ( v_norm2(v) >= MACHEPS*v_norm2(y)*A->m )
  361.     {
  362.     errmesg("spLUfactor()/spLUsolve()");
  363.     printf("# Sparse LU residual = %g [cf MACHEPS = %g]\n",
  364.            v_norm2(v), MACHEPS);
  365.     }
  366.     /* compute error in solution */
  367.     v_sub(x,z,z);
  368.     if ( v_norm2(z) > MACHEPS*v_norm2(x)*100*A->m )
  369.     {
  370.     errmesg("spLUfactor()/spLUsolve()");
  371.     printf("# Sparse LU solution error = %g [cf MACHEPS = %g]\n",
  372.            v_norm2(z), MACHEPS);
  373.     }
  374.  
  375.     /* now check spLUTsolve */
  376.     mem_stat_mark(4);
  377.     sp_vm_mlt(B,x,y);
  378.     spLUTsolve(A,pivot,y,z);
  379.     sp_vm_mlt(B,z,v);
  380.     mem_stat_free(4);
  381.  
  382.     /* compute residual */
  383.     v_sub(y,v,v);
  384.     if ( v_norm2(v) >= MACHEPS*v_norm2(y)*A->m )
  385.     {
  386.     errmesg("spLUTsolve()");
  387.     printf("# Sparse LU residual = %g [cf MACHEPS = %g]\n",
  388.            v_norm2(v), MACHEPS);
  389.     }
  390.     /* compute error in solution */
  391.     v_sub(x,z,z);
  392.     if ( v_norm2(z) > MACHEPS*v_norm2(x)*100*A->m )
  393.     {
  394.     errmesg("spLUTsolve()");
  395.     printf("# Sparse LU solution error = %g [cf MACHEPS = %g]\n",
  396.            v_norm2(z), MACHEPS);
  397.     }
  398.  
  399.     /* algebraic operations */
  400.     notice("addition,subtraction and multiplying by a number");
  401.     SP_FREE(A);
  402.     SP_FREE(B);
  403.  
  404.     m = 120;
  405.     n = 120;
  406.     deg = 5;
  407.     A = sp_get(m,n,deg);
  408.     B = sp_get(m,n,deg);
  409.     C = sp_get(m,n,deg);
  410.     C1 = sp_get(m,n,deg);
  411.  
  412.     for ( k = 0; k < m*deg; k++ )
  413.     {
  414.     i = (rand() >> 8) % m;
  415.     j = (rand() >> 8) % n;
  416.     sp_set_val(A,i,j,rand()/((Real)MAX_RAND));
  417.     i = (rand() >> 8) % m;
  418.     j = (rand() >> 8) % n;
  419.     sp_set_val(B,i,j,rand()/((Real)MAX_RAND));
  420.     }
  421.     
  422.     s1 = mrand(); 
  423.     B1 = sp_copy(B);
  424.  
  425.     mem_stat_mark(1);
  426.     sp_smlt(B,s1,C);
  427.     sp_add(A,C,C1);
  428.     sp_sub(C1,A,C);
  429.     sp_smlt(C,-1.0/s1,C);
  430.     sp_add(C,B1,C);
  431.  
  432.     s2 = 0.0;
  433.     for (k=0; k < C->m; k++) {
  434.        r = &(C->row[k]);
  435.        for (j=0; j < r->len; j++) {
  436.       if (s2 < fabs(r->elt[j].val)) 
  437.         s2 = fabs(r->elt[j].val);
  438.        }
  439.     }
  440.  
  441.     if (s2 > MACHEPS*A->m) {
  442.        errmesg("add, sub, mlt sparse matrices (args not in situ)\n");
  443.        printf(" difference = %g [MACEPS = %g]\n",s2,MACHEPS);
  444.     }
  445.  
  446.     sp_mltadd(A,B1,s1,C1);
  447.     sp_sub(C1,A,A);
  448.     sp_smlt(A,1.0/s1,C1);
  449.     sp_sub(C1,B1,C1);
  450.     mem_stat_free(1);
  451.  
  452.     s2 = 0.0;
  453.     for (k=0; k < C1->m; k++) {
  454.        r = &(C1->row[k]);
  455.        for (j=0; j < r->len; j++) {
  456.       if (s2 < fabs(r->elt[j].val)) 
  457.         s2 = fabs(r->elt[j].val);
  458.        }
  459.     }
  460.  
  461.     if (s2 > MACHEPS*A->m) {
  462.        errmesg("add, sub, mlt sparse matrices (args not in situ)\n");
  463.        printf(" difference = %g [MACEPS = %g]\n",s2,MACHEPS);
  464.     }
  465.  
  466.     V_FREE(x);
  467.     V_FREE(y);    
  468.     V_FREE(z);
  469.     V_FREE(u);
  470.     V_FREE(v);  
  471.     PX_FREE(pivot);
  472.     SP_FREE(A);
  473.     SP_FREE(B);
  474.     SP_FREE(C);
  475.     SP_FREE(B1);
  476.     SP_FREE(C1);
  477.  
  478.     printf("# Done testing (%s)\n",argv[0]);
  479.     mem_info();
  480. }
  481.     
  482.  
  483.  
  484.  
  485.  
  486.