home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / matop < prev    next >
Encoding:
Text File  |  1994-01-13  |  11.7 KB  |  526 lines

  1. initialises x so that x->ve[i] == i */
  2. VEC    *v_count(x)
  3. VEC    *x;
  4. {
  5.     int    i;
  6.  
  7.     if ( ! x )
  8.         error(E_NULL,"v_count");
  9.  
  10.     for ( i = 0; i < x->d printf("# det_mant1=%g, det_expt1=%d\n",
  11.         det_mant1,det_expt1); */
  12.      /* printf("# det_mant2=%g, det_expt2=%d\n",
  13.         det_mant2,det_expt2); */
  14.      if ( det_mant1 == 0.0 )
  15.      {   /* multiple e-val of T */
  16.         err_est->ve[i] = 0.0;
  17.         continue;
  18.      }
  19.      else if ( det_mant2 == 0.0 )
  20.      {
  21.         err_est->ve[i] = HUGE;
  22.         continue;
  23.      }
  24.      if ( (det_expt1 + det_expt2) % 2 )
  25.        /* if odd... */
  26.        det_mant = sqrt(2.0*fabs(det_mant1*det_mant2));
  27.      else /* if even... */
  28.        det_mant = sqrt(fabs(det_mant1*det_mant2));
  29.      det_expt = (det_expt1+det_expt2)/2;
  30.      err_est->ve[i] = fabs(beta*
  31.                    ldexp(pb_mant/det_mant,pb_expt-det_expt));
  32.       }
  33.    }
  34.    
  35.    return a;
  36. }
  37.  
  38. /* iter_splanczos2 -- version of iter_lanczos2() that uses sparse matrix data
  39.    structure */
  40.  
  41. VEC    *iter_splanczos2(A,m,x0,evals,err_est)
  42. SPMAT    *A;
  43. int     m;
  44. VEC    *x0;        /* initial vector */
  45. VEC    *evals;        /* eigenvalue vector */
  46. VEC    *err_est;    /* error estimates of eigenvalues */
  47. {    
  48.    ITER *ip;
  49.    VEC *a;
  50.    
  51.    ip = iter_get(0,0);
  52.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  53.    ip->A_par = (void *) A;
  54.    ip->x = x0;
  55.    ip->k = m;
  56.    a = iter_lanczos2(ip,evals,err_est);    
  57.    ip->shared_x = ip->shared_b = TRUE;
  58.    iter_free(ip);   /* release only ITER structure */
  59.    return a;
  60. }
  61.  
  62.  
  63.  
  64.  
  65. /*
  66.   Conjugate gradient method
  67.   Another variant - mainly for testing
  68.   */
  69.  
  70. VEC  *iter_cg1(ip)
  71. ITER *ip;
  72. {
  73.    static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
  74.    Real    alpha;
  75.    double inner,nres;
  76.    VEC *rr;   /* rr == r or rr == z */
  77.    
  78.    if (ip == INULL)
  79.      error(E_NULL,"iter_cg");
  80.    if (!ip->Ax || !ip->b)
  81.      error(E_NULL,"iter_cg");
  82.    if ( ip->x == ip->b )
  83.      error(E_INSITU,"iter_cg");
  84.    if (!ip->stop_crit)
  85.      error(E_NULL,"iter_cg");
  86.    
  87.    if ( ip->eps <= 0.0 )
  88.      ip->eps = MACHEPS;
  89.    
  90.    r = v_resize(r,ip->b->dim);
  91.    p = v_resize(p,ip->b->dim);
  92.    q = v_resize(q,ip->b->dim);
  93.    
  94.    MEM_STAT_REG(r,TYPE_VEC);
  95.    MEM_STAT_REG(p,TYPE_VEC);
  96.    MEM_STAT_REG(q,TYPE_VEC);
  97.    
  98.    if (ip->Bx != (Fun_Ax)NULL) {
  99.       z = v_resize(z,ip->b->dim);
  100.       MEM_STAT_REG(z,TYPE_VEC    error(E_SIZES,"m_sub");
  101.     if ( out==(MAT *)NULL || out->m != mat1->m || out->n != mat1->n )
  102.         out = m_resize(out,mat1->m,mat1->n);
  103.     m = mat1->m;    n = mat1->n;
  104.     for ( i=0; i<m; i++ )
  105.     {
  106.         __sub__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  107.         /**************************************************
  108.         for ( j=0; j<n; j++ )
  109.             out->me[i][j] = mat1->me[i][j]-mat2->me[i][j];
  110.         **************************************************/
  111.     }
  112.  
  113.     return (out);
  114. }
  115.  
  116. /* m_mlt -- matrix-matrix multiplication */
  117. MAT    *m_mlt(A,B,OUT)
  118. MAT    *A,*B,*OUT;
  119. {
  120.     u_int    i, /* j, */ k, m, n, p;
  121.     Real    **A_v, **B_v /*, *B_row, *OUT_row, sum, tmp */;
  122.  
  123.     if ( A==(MAT *)NULL || B==(MAT *)NULL )
  124.         error(E_NULL,"m_mlt");
  125.     if ( A->n != B->m )
  126.         error(E_SIZES,"m_mlt");
  127.     if ( A == OUT || B == OUT )
  128.         error(E_INSITU,"m_mlt");
  129.     m = A->m;    n = A->n;    p = B->n;
  130.     A_v = A->me;        B_v = B->me;
  131.  
  132.     if ( OUT==(MAT *)NULL || OUT->m != A->m || OUT->n != B->n )
  133.         OUT = m_resize(OUT,A->m,B->n);
  134.  
  135. /****************************************************************
  136.     for ( i=0; i<m; i++ )
  137.         for  ( j=0; j<p; j++ )
  138.         {
  139.             sum = 0.0;
  140.             for ( k=0; k<n; k++ )
  141.                 sum += A_v[i][k]*B_v[k][j];
  142.             OUT->me[i][j] = sum;
  143.         }
  144. ****************************************************************/
  145.     m_zero(OUT);
  146.     for ( i=0; i<m; i++ )
  147.         for ( k=0; k<n; k++ )
  148.         {
  149.             if ( A_v[i][k] != 0.0 )
  150.                 __mltadd__(OUT->me[i],B_v[k],A_v[i][k],(int)p);
  151.             /**************************************************
  152.             B_row = B_v[k];    OUT_row = OUT->me[i];
  153.             for ( j=0; j<p; j++ )
  154.             (*OUT_row++) += tmp*(*B_row++);
  155.             **************************************************/
  156.         }
  157.  
  158.     return OUT;
  159. }
  160.  
  161. /* mmtr_mlt -- matrix-matrix transposed multiplication
  162.     -- A.B^T is returned, and stored in OUT */
  163. MAT    *mmtr_mlt(A,B,OUT)
  164. MAT    *A, *B, *OUT;
  165. {
  166.     int    i, j, limit;
  167.     /* Real    *A_row, *B_row, sum; */
  168.  
  169.     if ( ! A || ! B )
  170.         error(E_NULL,"mmtr_mlt");
  171.     if ( A == OUT || B == OUT )
  172.         error(E_INSITU,"mmtr_mlt");
  173.     if ( A->n != B->n )
  174.         error(E_SIZES,"mmtr_mlt");
  175.     if ( ! OUT || OUT->m != A->m || OUT->n != B->m )
  176.         OUT = m_resize(OUT,A->m,B->m);
  177.  
  178.     limit = A->n;
  179.     for ( i = 0; i < A->m; i++ )
  180.         for ( j = 0; j < B->m; j++ )
  181.         {
  182.             OUT->me[i][j] = __ip__(A->me[i],B->me[j],(int)limit);
  183.             /**************************************************
  184.             sum = 0.0;
  185.             A_row = A->me[i];
  186.             B_row = B->me[j];
  187.             for ( k = 0; k < limit; k++ )
  188.             sum += (*A_row++)*(*B_row++);
  189.             OUT->me[i][j] = sum;
  190.             **************************************************/
  191.         }
  192.  
  193.     return OUT;
  194. }
  195.  
  196. /* mtrm_mlt -- matrix transposed-matrix multiplication
  197.     -- A^T.B is returned, result stored in OUT */
  198. MAT    *mtrm_mlt(A,B,OUT)
  199. MAT    *A, *B, *OUT;
  200. {
  201.     int    i, k, limit;
  202.     /* Real    *B_row, *OUT_row, multiplier; */
  203.  
  204.     if ( ! A || ! B )
  205.         error(E_NULL,"mmtr_mlt");
  206.     if ( A == OUT || B == OUT )
  207.         error(E_INSITU,"mtrm_mlt");
  208.     if ( A->m != B->m )
  209.         error(E_SIZES,"mmtr_mlt");
  210.     if ( ! OUT || OUT->m != A->n || OUT->n != B->n )
  211.         OUT = m_resize(OUT,A->n,B->n);
  212.  
  213.     limit = B->n;
  214.     m_zero(OUT);
  215.     for ( k = 0; k < A->m; k++ )
  216.         for ( i = 0; i < A->n; i++ )
  217.         {
  218.             if ( A->me[k][i] != 0.0 )
  219.             __mltadd__(OUT->me[i],B->me[k],A->me[k][i],(int)limit);
  220.             /**************************************************
  221.             multiplier = A->me[k][i];
  222.             OUT_row = OUT->me[i];
  223.             B_row   = B->me[k];
  224.             for ( j = 0; j < limit; j++ )
  225.             *(OUT_row++) += multiplier*(*B_row++);
  226.             **************************************************/
  227.         }
  228.  
  229.     return OUT;
  230. }
  231.  
  232. /* mv_mlt -- matrix-vector multiplication 
  233.         -- Note: b is treated as a column vector */
  234. VEC    *mv_mlt(A,b,out)
  235. MAT    *A;
  236. VEC    *b,*out;
  237. {
  238.     u_int    i, m, n;
  239.     Real    **A_v, *b_v /*, *A_row */;
  240.     /* register Real    sum; */
  241.  
  242.     if ( A==(MAT *)NULL || b==(VEC *)NULL )
  243.         error(E_NULL,"mv_mlt");
  244.     if ( A->n != b->dim )
  245.         error(E_SIZES,"mv_mlt");
  246.     if ( b == out )
  247.         error(E_INSITU,"mv_mlt");
  248.     if ( out == (VEC *)NULL || out->dim != A->m )
  249.         out = v_resize(out,A->m);
  250.  
  251.     m = A->m;        n = A->n;
  252.     A_v = A->me;        b_v = b->ve;
  253.     for ( i=0; i<m; i++ )
  254.     {
  255.         /* for ( j=0; j<n; j++ )
  256.             sum += A_v[i][j]*b_v[j]; */
  257.         out->ve[i] = __ip__(A_v[i],b_v,(int)n);
  258.         /**************************************************
  259.         A_row = A_v[i];        b_v = b->ve;
  260.         for ( j=0; j<n; j++ )
  261.             sum += (*A_row++)*(*b_v++);
  262.         out->ve[i] = sum;
  263.         **************************************************/
  264.     }
  265.  
  266.     return out;
  267. }
  268.  
  269. /* sm_mlt -- scalar-matrix multiply -- may be in-situ */
  270. MAT    *sm_mlt(scalar,matrix,out)
  271. double    scalar;
  272. MAT    *matrix,*out;
  273. {
  274.     u_int    m,n,i;
  275.  
  276.     if ( matrix==(MAT *)NULL )
  277.         error(E_NULL,"sm_mlt");
  278.     if ( out==(MAT *)NULL || out->m != matrix->m || out->n != matrix->n )
  279.         out = m_resize(out,matrix->m,matrix->n);
  280.     m = matrix->m;    n = matrix->n;
  281.     for ( i=0; i<m; i++ )
  282.         __smlt__(matrix->me[i],(double)scalar,out->me[i],(int)n);
  283.         /**************************************************
  284.         for ( j=0; j<n; j++ )
  285.             out->me[i][j] = scalar*matrix->me[i][j];
  286.         **************************************************/
  287.     return (out);
  288. }
  289.  
  290. /* vm_mlt -- vector-matrix multiplication 
  291.         -- Note: b is treated as a row vector */
  292. VEC    *vm_mlt(A,b,out)
  293. MAT    *A;
  294. VEC    *b,*out;
  295. {
  296.     u_int    j,m,n;
  297.     /* Real    sum,**A_v,*b_v; */
  298.  
  299.     if ( A==(MAT *)NULL || b==(VEC *)NULL )
  300.         error(E_NULL,"vm_mlt");
  301.     if ( A->m != b->dim )
  302.         error(E_SIZES,"vm_mlt");
  303.     if ( b == out )
  304.         error(E_INSITU,"vm_mlt");
  305.     if ( out == (VEC *)NULL || out->dim != A->n )
  306.         out = v_resize(out,A->n);
  307.  
  308.     m = A->m;        n = A->n;
  309.  
  310.     v_zero(out);
  311.     for ( j = 0; j < m; j++ )
  312.         if ( b->ve[j] != 0.0 )
  313.             __mltadd__(out->ve,A->me[j],b->ve[j],(int)n);
  314.     /**************************************************
  315.     A_v = A->me;        b_v = b->ve;
  316.     for ( j=0; j<n; j++ )
  317.     {
  318.         sum = 0.0;
  319.         for ( i=0; i<m; i++ )
  320.             sum += b_v[i]*A_v[i][j];
  321.         out->ve[j] = sum;
  322.     }
  323.     **************************************************/
  324.  
  325.     return out;
  326. }
  327.  
  328. /* m_transp -- transpose matrix */
  329. MAT    *m_transp(in,out)
  330. MAT    *in, *out;
  331. {
  332.     int    i, j;
  333.     int    in_situ;
  334.     Real    tmp;
  335.  
  336.     if ( in == (MAT *)NULL )
  337.         error(E_NULL,"m_transp");
  338.     if ( in == out && in->n != in->m )
  339.         error(E_INSITU2,"m_transp");
  340.     in_situ = ( in == out );
  341.     if ( out == (MAT *)NULL || out->m != in->n || out->n != in->m )
  342.         out = m_resize(out,in->n,in->m);
  343.  
  344.     if ( ! in_situ )
  345.         for ( i = 0; i < in->m; i++ )
  346.             for ( j = 0; j < in->n; j++ )
  347.                 out->me[j][i] = in->me[i][j];
  348.     else
  349.         for ( i = 1; i < in->m; i++ )
  350.             for ( j = 0; j < i; j++ )
  351.             {    tmp = in->me[i][j];
  352.                 in->me[i][j] = in->me[j][i];
  353.                 in->me[j][i] = tmp;
  354.             }
  355.  
  356.     return out;
  357. }
  358.  
  359. /* swap_rows -- swaps rows i and j of matrix A upto column lim */
  360. MAT    *swap_rows(A,i,j,lo,hi)
  361. MAT    *A;
  362. int    i, j, lo, hi;
  363. {
  364.     int    k;
  365.     Real    **A_me, tmp;
  366.  
  367.     if ( ! A )
  368.         error(E_NULL,"swap_rows");
  369.     if ( i < 0 || j < 0 || i >= A->m || j >= A->m )
  370.         error(E_SIZES,"swap_rows");
  371.     lo = max(0,lo);
  372.     hi = min(hi,A->n-1);
  373.     A_me = A->me;
  374.  
  375.     for ( k = lo; k <= hi; k++ )
  376.     {
  377.         tmp = A_me[k][i];
  378.         A_me[k][i] = A_me[k][j];
  379.         A_me[k][j] = tmp;
  380.     }
  381.     return A;
  382. }
  383.  
  384. /* swap_cols -- swap columns i and j of matrix A upto row lim */
  385. MAT    *swap_cols(A,i,j,lo,hi)
  386. MAT    *A;
  387. int    i, j, lo, hi;
  388. {
  389.     int    k;
  390.     Real    **A_me, tmp;
  391.  
  392.     if ( ! A )
  393.         error(E_NULL,"swap_cols");
  394.     if ( i < 0 || j < 0 || i >= A->n || j >= A->n )
  395.         error(E_SIZES,"swap_cols");
  396.     lo = max(0,lo);
  397.     hi = min(hi,A->m-1);
  398.     A_me = A->me;
  399.  
  400.     for ( k = lo; k <= hi; k++ )
  401.     {
  402.         tmp = A_me[i][k];
  403.         A_me[i][k] = A_me[j][k];
  404.         A_me[j][k] = tmp;
  405.     }
  406.     return A;
  407. }
  408.  
  409. /* ms_mltadd -- matrix-scalar multiply and add
  410.     -- may be in situ
  411.     -- returns out == A1 + s*A2 */
  412. MAT    *ms_mltadd(A1,A2,s,out)
  413. MAT    *A1, *A2, *out;
  414. double    s;
  415. {
  416.     /* register Real    *A1_e, *A2_e, *out_e; */
  417.     /* register int    j; */
  418.     int    i, m, n;
  419.  
  420.     if ( ! A1 || ! A2 )
  421.         error(E_NULL,"ms_mltadd");
  422.     if ( A1->m != A2->m || A1->n != A2->n )
  423.         error(E_SIZES,"ms_mltadd");
  424.  
  425.     if ( s == 0.0 )
  426.         return m_copy(A1,out);
  427.     if ( s == 1.0 )
  428.         return m_add(A1,A2,out);
  429.  
  430.     tracecatch(out = m_copy(A1,out),"ms_mltadd");
  431.  
  432.     m = A1->m;    n = A1->n;
  433.     for ( i = 0; i < m; i++ )
  434.     {
  435.         __mltadd__(out->me[i],A2->me[i],s,(int)n);
  436.         /**************************************************
  437.         A1_e = A1->me[i];
  438.         A2_e = A2->me[i];
  439.         out_e = out->me[i];
  440.         for ( j = 0; j < n; j++ )
  441.             out_e[j] = A1_e[j] + s*A2_e[j];
  442.         **************************************************/
  443.     }
  444.  
  445.     return out;
  446. }
  447.  
  448. /* mv_mltadd -- matrix-vector multiply and add
  449.     -- may not be in situ
  450.     -- returns out == v1 + alpha*A*v2 */
  451. VEC    *mv_mltadd(v1,v2,A,alpha,out)
  452. VEC    *v1, *v2, *out;
  453. MAT    *A;
  454. double    alpha;
  455. {
  456.     /* register    int    j; */
  457.     int    i, m, n;
  458.     Real    *v2_ve, *out_ve;
  459.  
  460.     if ( ! v1 || ! v2 || ! A )
  461.         error(E_NULL,"mv_mltadd");
  462.     if ( out == v2 )
  463.         error(E_INSITU,"mv_mltadd");
  464.     if ( v1->dim != A->m || v2->dim != A-> n )
  465.         error(E_SIZES,"mv_mltadd");
  466.  
  467.     tracecatch(out = v_copy(v1,out),"mv_mltadd");
  468.  
  469.     v2_ve = v2->ve;    out_ve = out->ve;
  470.     m = A->m;    n = A->n;
  471.  
  472.     if ( alpha == 0.0 )
  473.         return out;
  474.  
  475.     for ( i = 0; i < m; i++ )
  476.     {
  477.         out_ve[i] += alpha*__ip__(A->me[i],v2_ve,(int)n);
  478.         /**************************************************
  479.         A_e = A->me[i];
  480.         sum = 0.0;
  481.         for ( j = 0; j < n; j++ )
  482.             sum += A_e[j]*v2_ve[j];
  483.         out_ve[i] = v1->ve[i] + alpha*sum;
  484.         **************************************************/
  485.     }
  486.  
  487.     return out;
  488. }
  489.  
  490. /* vm_mltadd -- vector-matrix multiply and add
  491.     -- may not be in situ
  492.     -- returns out' == v1' + v2'*A */
  493. VEC    *vm_mltadd(v1,v2,A,alpha,out)
  494. VEC    *v1, *v2, *out;
  495. MAT    *A;
  496. double    alpha;
  497. {
  498.     int    /* i, */ j, m, n;
  499.     Real    tmp, /* *A_e, */ *out_ve;
  500.  
  501.     if ( ! v1 || ! v2 || ! A )
  502.         error(E_NULL,"vm_mltadd");
  503.     if ( v2 == out )
  504.         error(E_INSITU,"vm_mltadd");
  505.     if ( v1->dim != A->n || A->m != v2->dim )
  506.         error(E_SIZES,"vm_mltadd");
  507.  
  508.     tracecatch(out = v_copy(v1,out),"vm_mltadd");
  509.  
  510.     out_ve = out->ve;    m = A->m;    n = A->n;
  511.     for ( j = 0; j < m; j++ )
  512.     {
  513.         tmp = v2->ve[j]*alpha;
  514.         if ( tmp != 0.0 )
  515.             __mltadd__(out_ve,A->me[j],tmp,(int)n);
  516.         /**************************************************
  517.         A_e = A->me[j];
  518.         for ( i = 0; i < n; i++ )
  519.             out_ve[i] += A_e[i]*tmp;
  520.         **************************************************/
  521.     }
  522.  
  523.     return out;
  524. }
  525.  
  526.