home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / vecop < prev    next >
Encoding:
Text File  |  1994-03-08  |  13.1 KB  |  523 lines
  1.  [cf MACHEPS = %g]\n",
  2.            v_norm2(u), MACHEPS);
  3.     }
  4.  
  5.     Q = m_resize(Q,A->m,A->m);
  6.     makeQ(A,diag,Q);
  7.     makeR(A,A);
  8.     m_mlt(Q,A,C);
  9.     M_FREE(D);
  10.     D = m_get(B->m,B->n);
  11.     px_cols(pivot,C,D);
  12.     m_sub(B,D,D);
  13.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm1(B) )
  14.     {
  15.     errmesg("QRCPfactor()/makeQ()/makeR()");
  16.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  17.            m_norm1(D), MACHEPS);
  18.     }
  19.  
  20.     MEMCHK();
  21.  
  22.     /* Cholesky and LDL^T factorisation */
  23.     /* Use these for normal equations approach */
  24.     notice("Cholesky factor/solve");
  25.     mtrm_mlt(B,B,A);
  26.     CHfactor(A);
  27.     u = v_resize(u,B->n);
  28.     vm_mlt(B,y,u);
  29.     z = v_resize(z,B->n);
  30.     CHsolve(A,u,z);
  31.     v_sub(x,z,z);
  32.     if ( v_norm2(z) >= MACHEPS*v_norm2(x)*100 )
  33.     {
  34.     errmesg("CHfactor()/CHsolve()");
  35.     printf("# Cholesky solution error = %g [cf MACHEPS = %g]\n",
  36.            v_norm2(z), MACHEPS);
  37.     }
  38.     /* modified Cholesky factorisation should be identical with Cholesky
  39.        factorisation provided the matrix is "sufficiently positive definite */
  40.     mtrm_mlt(B,B,C);
  41.     MCHfactor(C,MACHEPS);
  42.     m_sub(A,C,C);
  43.     if ( m_norm1(C) >= MACHEPS*m_norm1(A) )
  44.     {
  45.     errmesg("MCHfactor()");
  46.     printf("# Modified Cholesky error = %g [cf MACHEPS = %g]\n",
  47.            m_norm1(C), MACHEPS);
  48.     }
  49.     /* now test the LDL^T factorisation -- using a negative def. matrix */
  50.     mtrm_mlt(B,B,A);
  51.     sm_mlt(-1.0,A,A);
  52.     m_copy(A,C);
  53.     LDLfactor(A);
  54.     LDLsolve(A,u,z);
  55.     w = v_get(A->m);
  56.     mv_mlt(C,z,w);
  57.     v_sub(w,u,w);
  58.     if ( v_norm2(w) >= MACHEPS*v_norm2(u)*m_norm1(C) )
  59.     {
  60.     errmesg("LDLfactor()/LDLsolve()");
  61.     printf("# LDL^T residual = %g [cf MACHEPS = %g]\n",
  62.            v_norm2(w), MACHEPS);
  63.     }
  64.     v_add(x,z,z);
  65.     if ( v_norm2(z) >= MACHEPS*v_norm2(x)*100 )
  66.     {
  67.     errmesg("LDLfactor()/LDLsolve()");
  68.     printf("# LDL^T solution error = %g [cf MACHEPS = %g]\n",
  69.            v_norm2(z), MACHEPS);
  70.     }
  71.  
  72.     MEMCHK();
  73.  
  74.     /* and now the Bunch-Kaufman-Parlett method */
  75.     /* set up D to be an indefinite diagonal matrix */
  76.     notice("Bunch-Kaufman-Parlett factor/solve");
  77.  
  78.     D = m_resize(D,B->m,B->m);
  79.     m_zero(D);
  80.     w = v_resize(w,B->m);
  81.     v_rand(w);
  82.     for ( i = 0; i < w->dim; i++ )
  83.     if ( v_entry(w,i) >= 0.5 )
  84.         m_set_val(D,i,i,1.0);
  85.     else
  86.         m_set_val(D,i,i,-1.0);
  87.     /* set A <- B^T.D.B */
  88.     C = m_resize(C,B->n,B->n);
  89.     C = mtrm_mlt(B,D,C);
  90.     A = m_mlt(C,B,A);
  91.     C = m_resize(C,B->n,B->n);
  92.     C = m_copy(A,C);
  93.     /* ... and use BKPfactor() */
  94.     blocks = px_get(A->m);
  95.     pivot = px_resize(pivot,A->m);
  96.     x = v_resize(x,A->m);
  97.     y = v_resize(y,A->m);
  98.     z = v_resize(z,A->m);
  99.     v_rand(x);
  100.     mv_mlt(A,x,y);
  101.     BKPfactor(A,pivot,blocks);
  102.     printf("# BKP pivot =\n");    px_output(pivot);
  103.     printf("# BKP blocks =\n");    px_output(blocks);
  104.     BKPsolve(A,pivot,blocks,y,z);
  105.     /* compute & check residual */
  106.     mv_mlt(C,z,w);
  107.     v_sub(w,y,w);
  108.     if ( v_norm2(w) >= MACHEPS*m_norm1(C)*v_norm2(z) )
  109.     {
  110.     errmesg("BKPfactor()/BKPsolve()");
  111.     printf("# BKP residual size = %g [cf MACHEPS = %g]\n",
  112.            v_norm2(w), MACHEPS);
  113.     }
  114.  
  115.     /* check update routines */
  116.     /* check LDLupdate() first */
  117.     notice("update L.D.L^T routine");
  118.     A = mtrm_mlt(B,B,A);
  119.     m_resize(C,A->m,A->n);
  120.     C = m_copy(A,C);
  121.     LDLfactor(A);
  122.     s1 = 3.7;
  123.     w = v_resize(w,A->m);
  124.     v_rand(w);
  125.     for ( i = 0; i < C->m; i++ )
  126.     for ( j = 0; j < C->n; j++ )
  127.         m_set_val(C,i,j,m_entry(C,i,j)+s1*v_entry(w,i)*v_entry(w,j));
  128.     LDLfactor(C);
  129.     LDLupdate(A,w,s1);
  130.     /* zero out strictly upper triangular parts of A and C */
  131.     for ( i = 0; i < A->m; i++ )
  132.     for ( j = i+1; j < A->n; j++ )
  133.     {
  134.         m_set_val(A,i,j,0.0);
  135.         m_set_val(C,i,j,0.0);
  136.     }
  137.     if ( m_norm1(m_sub(A,C,C)) >= sqrt(MACHEPS)*m_norm1(A) )
  138.     {
  139.     errmesg("LDLupdate()");
  140.     printf("# LDL update matrix error = %g [cf MACHEPS = %g]\n",
  141.            m_norm1(C), MACHEPS);
  142.     }
  143.  
  144.  
  145.     /* BAND MATRICES */
  146.  
  147. #define COL 40
  148. #define UDIAG  5
  149. #define LDIAG  2
  150.  
  151.    smrand(101);
  152.    bA = bd_get(LDIAG,UDIAG,COL);
  153.    bB = bd_get(LDIAG,UDIAG,COL);
  154.    bC = bd_get(LDIAG,UDIAG,COL);
  155.    A = m_resize(A,COL,COL);
  156.    B = m_resize(B,COL,COL);
  157.    pivot = px_resize(pivot,COL);
  158.    x = v_resize(x,COL);
  159.    w = v_resize(w,COL);
  160.    z = v_resize(z,COL);
  161.  
  162.    m_rand(A); 
  163.    /* generate band matrix */
  164.    mat2band(A,LDIAG,UDIAG,bA);
  165.    band2mat(bA,A);    /* now A is banded */
  166.    bB = bd_copy(bA,bB); 
  167.  
  168.    v_rand(x);  
  169.    mv_mlt(A,x,w);
  170.    z = v_copy(w,z);
  171.  
  172.    notice("band LU factorization");
  173.    bdLUfactor(bA,pivot);
  174.  
  175.    /* pivot will be changed */
  176.    bdLUsolve(bA,pivot,z,z);
  177.    v_sub(x,z,z);
  178.    if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
  179.       errmesg("incorrect solution (band LU factorization)");
  180.       printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
  181.          v_norm2(z),MACHEPS);
  182.    }
  183.  
  184.    /* solve transpose system */
  185.  
  186.    notice("band LU factorization for transpose system");
  187.    m_transp(A,B);
  188.    mv_mlt(B,x,w);
  189.  
  190.    bd_copy(bB,bA);
  191.    bd_transp(bA,bA);  
  192.    /* transposition in situ */
  193.    bd_transp(bA,bA);
  194.    bd_transp(bA,bB);
  195.  
  196.    bdLUfactor(bB,pivot);
  197.  
  198.    bdLUsolve(bB,pivot,w,z);
  199.    v_sub(x,z,z);
  200.    if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
  201.       errmesg("incorrect solution (band transposed LU factorization)");
  202.       printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
  203.          v_norm2(z),MACHEPS);
  204.    }
  205.  
  206.  
  207.    /* Cholesky factorization */
  208.  
  209.    notice("band Choleski LDL' factorization");
  210.    m_add(A,B,A);  /* symmetric matrix */
  211.    for (i=0; i < COL; i++)     /* positive definite */
  212.      A->me[i][i] += 2*LDIAG;   
  213.  
  214.    mat2band(A,LDIAG,LDIAG,bA);
  215.    band2mat(bA,A);              /* corresponding matrix A */
  216.  
  217.    v_rand(x);
  218.    mv_mlt(A,x,w);
  219.    z = v_copy(w,z);
  220.    
  221.    bdLDLfactor(bA);
  222.  
  223.    z = bdLDLsolve(bA,z,z);
  224.    v_sub(x,z,z);
  225.    if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
  226.       errmesg("incorrect solution (band LDL' factorization)");
  227.       printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
  228.          v_norm2(z),MACHEPS);
  229.    }
  230.  
  231.    /* new bandwidths */
  232.    m_rand(A);
  233.    bA = bd_resize(bA,UDIAG,LDIAG,COL);
  234.    bB = bd_resize(bB,UDIAG,LDIAG,COL);
  235.    mat2band(A,UDIAG,LDIAG,bA);
  236.    band2mat(bA,A);
  237.    bd_copy(bA,bB);
  238.  
  239.    mv_mlt(A,x,w);
  240.  
  241.    notice("band LU factorization (resized)");
  242.    bdLUfactor(bA,pivot);
  243.  
  244.    /* pivot will be changed */
  245.    bdLUsolve(bA,pivot,w,z);
  246.    v_sub(x,z,z);
  247.    if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
  248.       errmesg("incorrect solution (band LU factorization)");
  249.       printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
  250.          v_norm2(z),MACHEPS);
  251.    }
  252.  
  253.    /* testing transposition */
  254.  
  255.    notice("band matrix transposition");
  256.    m_zero(bA->mat);
  257.    bd_copy(bB,bA);
  258.    m_zero(bB->mat);
  259.    bd_copy(bA,bB);
  260.  
  261.    bd_transp(bB,bB);
  262.    bd_transp(bB,bB);
  263.  
  264.    m_zero(bC->mat);
  265.    bd_copy(bB,bC);
  266.  
  267.    m_sub(bA->mat,bC->mat,bC->mat);
  268.    if (m_norm_inf(bC->mat) > MACHEPS*bC->mat->n) {
  269.       errmesg("band transposition");
  270.       printf(" difference ||A - (A')'|| = %g\n",m_norm_inf(bC->mat));
  271.    }
  272.  
  273.    bd_free(bA);
  274.    bd_free(bB);
  275.    bd_free(bC);
  276.  
  277.  
  278.     MEMCHK();
  279.  
  280.     /* now check QRupdate() routine */
  281.     notice("update QR routine");
  282.  
  283.     B = m_resize(B,15,7);
  284.     A = m_resize(A,B->m,B->n);
  285.     m_copy(B,A);
  286.     diag = v_resize(diag,A->n);
  287.     beta = v_resize(beta,A->n);
  288.     QRfactor(A,diag);
  289.     Q = m_resize(Q,A->m,A->m);
  290.     makeQ(A,diag,Q);
  291.     makeR(A,A);
  292.     m_resize(C,A->m,A->n);
  293.     w = v_resize(w,A->m);
  294.     v = v_resize(v,A->n);
  295.     u = v_resize(u,A->m);
  296.     v_rand(w);
  297.     v_rand(v);
  298.     vm_mlt(Q,w,u);
  299.     QRupdate(Q,A,u,v);
  300.     m_mlt(Q,A,C);
  301.     for ( i = 0; i < B->m; i++ )
  302.     for ( j = 0; j < B->n; j++ )
  303.         m_set_val(B,i,j,m_entry(B,i,j)+v_entry(w,i)*v_entry(v,j));
  304.     m_sub(B,C,C);
  305.     if ( m_norm1(C) >= MACHEPS*m_norm1(A)*m_norm1(Q)*2 )
  306.     {
  307.     errmesg("QRupdate()");
  308.     printf("# Reconstruction error in QR update = %g [cf MACHEPS = %g]\n",
  309.            m_norm1(C), MACHEPS);
  310.     }
  311.     m_resize(D,Q->m,Q->n);
  312.     mtrm_mlt(Q,Q,D);
  313.     for ( i = 0; i < D->m; i++ )
  314.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  315.     if ( m_norm1(D) >= 10*MACHEPS*m_norm1(Q)*m_norm_inf(Q) )
  316.     {
  317.     errmesg("QRupdate()");
  318.     printf("# QR update orthogonality error = %g [cf MACHEPS = %g]\n",
  319.            m_norm1(D), MACHEPS);
  320.     }
  321.  
  322.     /* Now check eigenvalue/SVD routines */
  323.     notice("eigenvalue and SVD routines");
  324.     A = m_resize(A,11,11);
  325.     B = m_resize(B,A->m,A->n);
  326.     C = m_resize(C,A->m,A->n);
  327.     D = m_resize(D,A->m,A->n);
  328.     Q = m_resize(Q,A->m,A->n);
  329.  
  330.     m_rand(A);
  331.     /* A <- A + A^T  for symmetric case */
  332.     m_add(A,m_transp(A,C),A);
  333.     u = v_resize(u,A->m);
  334.     u = symmeig(A,Q,u);
  335.     m_zero(B);
  336.     for ( i = 0; i < B->m; i++ )
  337.     m_set_val(B,i,i,v_entry(u,i));
  338.     m_mlt(Q,B,C);
  339.     mmtr_mlt(C,Q,D);
  340.     m_sub(A,D,D);
  341.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*v_norm_inf(u)*3 )
  342.     {
  343.     errmesg("symmeig()");
  344.     printf("# Reconstruction error = %g [cf MACHEPS = %g]\n",
  345.            m_norm1(D), MACHEPS);
  346.     }
  347.     mtrm_mlt(Q,Q,D);
  348.     for ( i = 0; i < D->m; i++ )
  349.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  350.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*3 )
  351.     {
  352.     errmesg("symmeig()");
  353.     printf("# symmeig() orthogonality error = %g [cf MACHEPS = %g]\n",
  354.            m_norm1(D), MACHEPS);
  355.     }
  356.  
  357.     MEMCHK();
  358.  
  359.     /* now test (real) Schur decomposition */
  360.     /* m_copy(A,B); */
  361.     M_FREE(A);
  362.     A = m_get(11,11);
  363.     m_rand(A);
  364.     B = m_copy(A,B);
  365.     MEMCHK();
  366.  
  367.     B = schur(B,Q);
  368.     MEMCHK();
  369.  
  370.     m_mlt(Q,B,C);
  371.     mmtr_mlt(C,Q,D);
  372.     MEMCHK();
  373.     m_sub(A,D,D);
  374.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*m_norm1(B)*5 )
  375.     {
  376.     errmesg("schur()");
  377.     printf("# Schur reconstruction error = %g [cf MACHEPS = %g]\n",
  378.            m_norm1(D), MACHEPS);
  379.     }
  380.  
  381.     /* orthogonality check */
  382.     mmtr_mlt(Q,Q,D);
  383.     for ( i = 0; i < D->m; i++ )
  384.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  385.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*10 )
  386.     {
  387.     errmesg("schur()");
  388.     printf("# Schur orthogonality error = %g [cf MACHEPS = %g]\n",
  389.            m_norm1(D), MACHEPS);
  390.     }
  391.  
  392.     MEMCHK();
  393.  
  394.     /* check for termination */
  395.     Atmp = m_get(2,2);
  396.     Qtmp = m_get(2,2);
  397.     /* this is a 2 x 2 matrix with real eigenvalues */
  398.     Atmp->me[0][0] = 1;
  399.     Atmp->me[1][1] = 1;
  400.     Atmp->me[0][1] = 4;
  401.     Atmp->me[1][0] = 1;
  402.     schur(Atmp,Qtmp);
  403.  
  404.     MEMCHK();
  405.  
  406.     /* now test SVD */
  407.     A = m_resize(A,11,7);
  408.     m_rand(A);
  409.     U = m_get(A->n,A->n);
  410.     Q = m_resize(Q,A->m,A->m);
  411.     u = v_resize(u,max(A->m,A->n));
  412.     svd(A,Q,U,u);
  413.     /* check reconstruction of A */
  414.     D = m_resize(D,A->m,A->n);
  415.     C = m_resize(C,A->m,A->n);
  416.     m_zero(D);
  417.     for ( i = 0; i < min(A->m,A->n); i++ )
  418.     m_set_val(D,i,i,v_entry(u,i));
  419.     mtrm_mlt(Q,D,C);
  420.     m_mlt(C,U,D);
  421.     m_sub(A,D,D);
  422.     if ( m_norm1(D) >= MACHEPS*m_norm1(U)*m_norm_inf(Q)*m_norm1(A) )
  423.     {
  424.     errmesg("svd()");
  425.     printf("# SVD reconstruction error = %g [cf MACHEPS = %g]\n",
  426.            m_norm1(D), MACHEPS);
  427.     }
  428.     /* check orthogonality of Q and U */
  429.     D = m_resize(D,Q->n,Q->n);
  430.     mtrm_mlt(Q,Q,D);
  431.     for ( i = 0; i < D->m; i++ )
  432.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  433.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*5 )
  434.     {
  435.     errmesg("svd()");
  436.     printf("# SVD orthognality error (Q) = %g [cf MACHEPS = %g\n",
  437.            m_norm1(D), MACHEPS);
  438.     }
  439.     D = m_resize(D,U->n,U->n);
  440.     mtrm_mlt(U,U,D);
  441.     for ( i = 0; i < D->m; i++ )
  442.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  443.     if ( m_norm1(D) >= MACHEPS*m_norm1(U)*m_norm_inf(U)*5 )
  444.     {
  445.     errmesg("svd()");
  446.     printf("# SVD orthognality error (U) = %g [cf MACHEPS = %g\n",
  447.            m_norm1(D), MACHEPS);
  448.     }
  449.     for ( i = 0; i < u->dim; i++ )
  450.     if ( v_entry(u,i) < 0 || (i < u->dim-1 &&
  451.                   v_entry(u,i+1) > v_entry(u,i)) )
  452.         break;
  453.     if ( i < u->dim )
  454.     {
  455.     errmesg("svd()");
  456.     printf("# SVD sorting error\n");
  457.     }
  458.  
  459.  
  460.     /* test of long vectors */
  461.     notice("Long vectors");
  462.     x = v_resize(x,100000);
  463.     y = v_resize(y,100000);
  464.     z = v_resize(z,100000);
  465.     v_rand(x);
  466.     v_rand(y);
  467.     v_mltadd(x,y,3.0,z);
  468.     sv_mlt(1.0/3.0,z,z);
  469.     v_mltadd(z,x,-1.0/3.0,z);
  470.     v_sub(z,y,x);
  471.     if (v_norm2(x) >= MACHEPS*(x->dim)) {
  472.        errmesg("long vectors");
  473.        printf(" norm = %g\n",v_norm2(x));
  474.     }
  475.  
  476.     mem_stat_free(1);
  477.  
  478.     MEMCHK();
  479.  
  480.     /**************************************************
  481.     VEC        *x, *y, *z, *u, *v, *w;
  482.     VEC        *diag, *beta;
  483.     PERM    *pi1, *pi2, *pi3, *pivot, *blocks;
  484.     MAT        *A, *B, *C, *D, *Q, *U;
  485.     **************************************************/
  486.     V_FREE(x);        V_FREE(y);    V_FREE(z);
  487.     V_FREE(u);        V_FREE(v);    V_FREE(w);
  488.     V_FREE(diag);    V_FREE(beta);
  489.     PX_FREE(pi1);    PX_FREE(pi2);    PX_FREE(pi3);
  490.     PX_FREE(pivot);    PX_FREE(blocks);
  491.     M_FREE(A);        M_FREE(B);    M_FREE(C);
  492.     M_FREE(D);        M_FREE(Q);    M_FREE(U);
  493.     M_FREE(Atmp);    M_FREE(Qtmp);
  494.  
  495.     MEMCHK();
  496.     printf("# Finished torture test\n");
  497.     mem_info();
  498.  
  499.     return 0;
  500. }
  501.  
  502. FileDataŵtutadv†Eÿÿÿ@ò XÂ
  503. /* routines from the section 8 of tutorial.txt */
  504.  
  505. #include "matrix.h"
  506.  
  507. #define M3D_LIST    3      /* list number */
  508. #define TYPE_MAT3D  0      /* the number of a type */
  509.  
  510. /* type for 3 dimensional matrices */
  511. typedef struct {
  512.     int l,m,n;    /* actual dimensions */
  513.     int max_l, max_m, max_n;    /* maximal dimensions */
  514.     Real ***me;    /* pointer to matrix elements */
  515.                    /* we do not consider segmented memory */
  516.         Real *base, **me2d;  /* me and me2d are additional pointers 
  517.                 to base */
  518. } MAT3D;
  519.  
  520.  
  521. /* function for creating a variable of MAT3D type */
  522.  
  523. MAT3D