home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / memstat < prev    next >
Encoding:
Text File  |  1994-01-13  |  8.4 KB  |  366 lines
  1. ,B,OUT)
  2. MAT    *A,*B,*OUT;
  3. {
  4.     u_int    i, /* j, */ k, m, n, p;
  5.     Real    **A_v, **B_v /*, *B_row, *OUT_row, sum, tmp */;
  6.  
  7.     if ( A==(MAT *)NULL || B==(MAT *)NULL )
  8.         error(E_NULL,"m_mlt");
  9.     if ( A->n != B->m )
  10.         error(E_SIZES,"m_mlt");
  11.     if ( A == OUT || B == OUT )
  12.         error(E_INSITU,"m_mlt");
  13.     m = A->m;    n = A->n;    p = B->n;
  14.     A_v = A->me;        B_v = B->me;
  15.  
  16.     if ( OUT==(MAT *)NULL || OUT->m != A->m || OUT->n != B->n )
  17.         OUT = m_resize(OUT,A->m,B->n);
  18.  
  19. /****************************************************************
  20.     for ( i=0; i<m; i++ )
  21.         for  ( j=0; j<p; j++ )
  22.         {
  23.             sum = 0.0;
  24.             for ( k=0; k<n; k++ )
  25.                 sum += A_v[i][k]*B_v[k][j];
  26.             OUT->me[i][j] = sum;
  27.         }
  28. ****************************************************************/
  29.     m_zero(OUT);
  30.     for ( i=0; i<m; i++ )
  31.         for ( k=0; k<n; k++ )
  32.         {
  33.             if ( A_v[i][k] != 0.0 )
  34.                 __mltadd__(OUT->me[i],B_v[k],A_v[i][k],(int)p);
  35.             /**************************************************
  36.             B_row = B_v[k];    OUT_row = OUT->me[i];
  37.             for ( j=0; j<p; j++ )
  38.             (*OUT_row++) += tmp*(*B_row++);
  39.             **************************************************/
  40.         }
  41.  
  42.     return OUT;
  43. }
  44.  
  45. /* mmtr_mlt -- matrix-matrix transposed multiplication
  46.     -- A.B^T is returned, and stored in OUT */
  47. MAT    *mmtr_mlt(A,B,OUT)
  48. MAT    *A, *B, *OUT;
  49. {
  50.     int    i, j, limit;
  51.     /* Real    *A_row, *B_row, sum; */
  52.  
  53.     if ( ! A || ! B )
  54.         error(E_NULL,"mmtr_mlt");
  55.     if ( A == OUT || B == OUT )
  56.         error(E_INSITU,"mmtr_mlt");
  57.     if ( A->n != B->n )
  58.         error(E_SIZES,"mmtr_mlt");
  59.     if ( ! OUT || OUT->m != A->m || OUT->n != B->m )
  60.         OUT = m_resize(OUT,A->m,B->m);
  61.  
  62.     limit = A->n;
  63.     for ( i = 0; i < A->m; i++ )
  64.         for ( j = 0; j < B->m; j++ )
  65.         {
  66.             OUT->me[i][j] = __ip__(A->me[i],B->me[j],(int)limit);
  67.             /**************************************************
  68.             sum = 0.0;
  69.             A_row = A->me[i];
  70.             B_row = B->me[j];
  71.             for ( k = 0; k < limit; k++ )
  72.             sum += (*A_row++)*(*B_row++);
  73.             OUT->me[i][j] = sum;
  74.             **************************************************/
  75.         }
  76.  
  77.     return OUT;
  78. }
  79.  
  80. /* mtrm_mlt -- matrix transposed-matrix multiplication
  81.     -- A^T.B is returned, result stored in OUT */
  82. MAT    *mtrm_mlt(A,B,OUT)
  83. MAT    *A, *B, *OUT;
  84. {
  85.     int    i, k, limit;
  86.     /* Real    *B_row, *OUT_row, multiplier; */
  87.  
  88.     if ( ! A || ! B )
  89.         error(E_NULL,"mmtr_mlt");
  90.     if ( A == OUT || B == OUT )
  91.         error(E_INSITU,"mtrm_mlt");
  92.     if ( A->m != B->m )
  93.         error(E_SIZES,"mmtr_mlt");
  94.     if ( ! OUT || OUT->m != A->n || OUT->n != B->n )
  95.         OUT = m_resize(OUT,A->n,B->n);
  96.  
  97.     limit = B->n;
  98.     m_zero(OUT);
  99.     for ( k = 0; k < A->m; k++ )
  100.         for ( i = 0; i < A->n; i++ )
  101.         {
  102.             if ( A->me[k][i] != 0.0 )
  103.             __mltadd__(OUT->me[i],B->me[k],A->me[k][i],(int)limit);
  104.             /**************************************************
  105.             multiplier = A->me[k][i];
  106.             OUT_row = OUT->me[i];
  107.             B_row   = B->me[k];
  108.             for ( j = 0; j < limit; j++ )
  109.             *(OUT_row++) += multiplier*(*B_row++);
  110.             **************************************************/
  111.         }
  112.  
  113.     return OUT;
  114. }
  115.  
  116. /* mv_mlt -- matrix-vector multiplication 
  117.         -- Note: b is treated as a column vector */
  118. VEC    *mv_mlt(A,b,out)
  119. MAT    *A;
  120. VEC    *b,*out;
  121. {
  122.     u_int    i, m, n;
  123.     Real    **A_v, *b_v /*, *A_row */;
  124.     /* register Real    sum; */
  125.  
  126.     if ( A==(MAT *)NULL || b==(VEC *)NULL )
  127.         error(E_NULL,"mv_mlt");
  128.     if ( A->n != b->dim )
  129.         error(E_SIZES,"mv_mlt");
  130.     if ( b == out )
  131.         error(E_INSITU,"mv_mlt");
  132.     if ( out == (VEC *)NULL || out->dim != A->m )
  133.         out = v_resize(out,A->m);
  134.  
  135.     m = A->m;        n = A->n;
  136.     A_v = A->me;        b_v = b->ve;
  137.     for ( i=0; i<m; i++ )
  138.     {
  139.         /* for ( j=0; j<n; j++ )
  140.             sum += A_v[i][j]*b_v[j]; */
  141.         out->ve[i] = __ip__(A_v[i],b_v,(int)n);
  142.         /**************************************************
  143.         A_row = A_v[i];        b_v = b->ve;
  144.         for ( j=0; j<n; j++ )
  145.             sum += (*A_row++)*(*b_v++);
  146.         out->ve[i] = sum;
  147.         **************************************************/
  148.     }
  149.  
  150.     return out;
  151. }
  152.  
  153. /* sm_mlt -- scalar-matrix multiply -- may be in-situ */
  154. MAT    *sm_mlt(scalar,matrix,out)
  155. double    scalar;
  156. MAT    *matrix,*out;
  157. {
  158.     u_int    m,n,i;
  159.  
  160.     if ( matrix==(MAT *)NULL )
  161.         error(E_NULL,"sm_mlt");
  162.     if ( out==(MAT *)NULL || out->m != matrix->m || out->n != matrix->n )
  163.         out = m_resize(out,matrix->m,matrix->n);
  164.     m = matrix->m;    n = matrix->n;
  165.     for ( i=0; i<m; i++ )
  166.         __smlt__(matrix->me[i],(double)scalar,out->me[i],(int)n);
  167.         /**************************************************
  168.         for ( j=0; j<n; j++ )
  169.             out->me[i][j] = scalar*matrix->me[i][j];
  170.         **************************************************/
  171.     return (out);
  172. }
  173.  
  174. /* vm_mlt -- vector-matrix multiplication 
  175.         -- Note: b is treated as a row vector */
  176. VEC    *vm_mlt(A,b,out)
  177. MAT    *A;
  178. VEC    *b,*out;
  179. {
  180.     u_int    j,m,n;
  181.     /* Real    sum,**A_v,*b_v; */
  182.  
  183.     if ( A==(MAT *)NULL || b==(VEC *)NULL )
  184.         error(E_NULL,"vm_mlt");
  185.     if ( A->m != b->dim )
  186.         error(E_SIZES,"vm_mlt");
  187.     if ( b == out )
  188.         error(E_INSITU,"vm_mlt");
  189.     if ( out == (VEC *)NULL || out->dim != A->n )
  190.         out = v_resize(out,A->n);
  191.  
  192.     m = A->m;        n = A->n;
  193.  
  194.     v_zero(out);
  195.     for ( j = 0; j < m; j++ )
  196.         if ( b->ve[j] != 0.0 )
  197.             __mltadd__(out->ve,A->me[j],b->ve[j],(int)n);
  198.     /**************************************************
  199.     A_v = A->me;        b_v = b->ve;
  200.     for ( j=0; j<n; j++ )
  201.     {
  202.         sum = 0.0;
  203.         for ( i=0; i<m; i++ )
  204.             sum += b_v[i]*A_v[i][j];
  205.         out->ve[j] = sum;
  206.     }
  207.     **************************************************/
  208.  
  209.     return out;
  210. }
  211.  
  212. /* m_transp -- transpose matrix */
  213. MAT    *m_transp(in,out)
  214. MAT    *in, *out;
  215. {
  216.     int    i, j;
  217.     int    in_situ;
  218.     Real    tmp;
  219.  
  220.     if ( in == (MAT *)NULL )
  221.         error(E_NULL,"m_transp");
  222.     if ( in == out && in->n != in->m )
  223.         error(E_INSITU2,"m_transp");
  224.     in_situ = ( in == out );
  225.     if ( out == (MAT *)NULL || out->m != in->n || out->n != in->m )
  226.         out = m_resize(out,in->n,in->m);
  227.  
  228.     if ( ! in_situ )
  229.         for ( i = 0; i < in->m; i++ )
  230.             for ( j = 0; j < in->n; j++ )
  231.                 out->me[j][i] = in->me[i][j];
  232.     else
  233.         for ( i = 1; i < in->m; i++ )
  234.             for ( j = 0; j < i; j++ )
  235.             {    tmp = in->me[i][j];
  236.                 in->me[i][j] = in->me[j][i];
  237.                 in->me[j][i] = tmp;
  238.             }
  239.  
  240.     return out;
  241. }
  242.  
  243. /* swap_rows -- swaps rows i and j of matrix A upto column lim */
  244. MAT    *swap_rows(A,i,j,lo,hi)
  245. MAT    *A;
  246. int    i, j, lo, hi;
  247. {
  248.     int    k;
  249.     Real    **A_me, tmp;
  250.  
  251.     if ( ! A )
  252.         error(E_NULL,"swap_rows");
  253.     if ( i < 0 || j < 0 || i >= A->m || j >= A->m )
  254.         error(E_SIZES,"swap_rows");
  255.     lo = max(0,lo);
  256.     hi = min(hi,A->n-1);
  257.     A_me = A->me;
  258.  
  259.     for ( k = lo; k <= hi; k++ )
  260.     {
  261.         tmp = A_me[k][i];
  262.         A_me[k][i] = A_me[k][j];
  263.         A_me[k][j] = tmp;
  264.     }
  265.     return A;
  266. }
  267.  
  268. /* swap_cols -- swap columns i and j of matrix A upto row lim */
  269. MAT    *swap_cols(A,i,j,lo,hi)
  270. MAT    *A;
  271. int    i, j, lo, hi;
  272. {
  273.     int    k;
  274.     Real    **A_me, tmp;
  275.  
  276.     if ( ! A )
  277.         error(E_NULL,"swap_cols");
  278.     if ( i < 0 || j < 0 || i >= A->n || j >= A->n )
  279.         error(E_SIZES,"swap_cols");
  280.     lo = max(0,lo);
  281.     hi = min(hi,A->m-1);
  282.     A_me = A->me;
  283.  
  284.     for ( k = lo; k <= hi; k++ )
  285.     {
  286.         tmp = A_me[i][k];
  287.         A_me[i][k] = A_me[j][k];
  288.         A_me[j][k] = tmp;
  289.     }
  290.     return A;
  291. }
  292.  
  293. /* ms_mltadd -- matrix-scalar multiply and add
  294.     -- may be in situ
  295.     -- returns out == A1 + s*A2 */
  296. MAT    *ms_mltadd(A1,A2,s,out)
  297. MAT    *A1, *A2, *out;
  298. double    s;
  299. {
  300.     /* register Real    *A1_e, *A2_e, *out_e; */
  301.     /* register int    j; */
  302.     int    i, m, n;
  303.  
  304.     if ( ! A1 || ! A2 )
  305.         error(E_NULL,"ms_mltadd");
  306.     if ( A1->m != A2->m || A1->n != A2->n )
  307.         error(E_SIZES,"ms_mltadd");
  308.  
  309.     if ( s == 0.0 )
  310.         return m_copy(A1,out);
  311.     if ( s == 1.0 )
  312.         return m_add(A1,A2,ou2(v_sub(x,y,z)) >= MACHEPS )
  313.      errmesg("PERMute vector");
  314.    pi2 = px_inv(pi1,pi2);
  315.    pi3 = px_mlt(pi1,pi2,pi3);
  316.    for ( i = 0; i < pi3->size; i++ )
  317.      if ( pi3->pe[i] != i )
  318.        errmesg("PERM inverse/multiply");
  319.    
  320.    n = px_free_vars(&pi1,&pi2,&pi3,NULL);
  321.    if ( n != 3 )
  322.      errmesg("PERM px_free_vars"); 
  323.  
  324. #ifdef SPARSE   
  325.    /* set up two random sparse matrices */
  326.    m = 120;
  327.    n = 100;
  328.    deg = 5;
  329.    notice("allocating sparse matrices");
  330.    k = sp_get_vars(m,n,deg,&sA,&sB,NULL);
  331.    if (k != 2) {
  332.       errmesg("sp_get_vars");
  333.       printf(" n = %d (should be 2)\n",k);
  334.    }
  335.    
  336.    notice("setting and getting matrix entries");
  337.    for ( k = 0; k < m*deg; k++ )
  338.    {
  339.       i = (rand() >> 8) % m;
  340.       j = (rand() >> 8) % n;
  341.       sp_set_val(sA,i,j,rand()/((Real)MAX_RAND));
  342.       i = (rand() >> 8) % m;
  343.       j = (rand() >> 8) % n;
  344.       sp_set_val(sB,i,j,rand()/((Real)MAX_RAND));
  345.    }
  346.    for ( k = 0; k < 10; k++ )
  347.    {
  348.       s1 = rand()/((Real)MAX_RAND);
  349.       i = (rand() >> 8) % m;
  350.       j = (rand() >> 8) % n;
  351.       sp_set_val(sA,i,j,s1);
  352.       s2 = sp_get_val(sA,i,j);
  353.       if ( fabs(s1 - s2) >= MACHEPS ) {
  354.      printf(" s1 = %g, s2 = %g, |s1 - s2| = %g\n", 
  355.         s1,s2,fabs(s1-s2));
  356.      break;
  357.       }
  358.    }
  359.    if ( k < 10 )
  360.      errmesg("sp_set_val()/sp_get_val()");
  361.    
  362.    /* check column access paths */
  363.    notice("resizing and access paths");
  364.    k = sp_resize_vars(sA->m+10,sA->n+10,&sA,&sB,NULL);
  365.    if (k != 2) {
  366.       errmesg("sp_ge