home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / DOC / fnindex next >
Encoding:
Text File  |  1994-01-13  |  16.8 KB  |  604 lines
  1. );
  2.  
  3.     mem_stat_mark(1);
  4.  
  5.     notice("complex arithmetic & special functions");
  6.  
  7.     ONE = zmake(1.0,0.0);
  8.     printf("# ONE = ");    z_output(ONE);
  9.     z1.re = mrand();    z1.im = mrand();
  10.     z2.re = mrand();    z2.im = mrand();
  11.     z3 = zadd(z1,z2);
  12.     if ( fabs(z1.re+z2.re-z3.re) + fabs(z1.im+z2.im-z3.im) > 10*MACHEPS )
  13.     errmesg("zadd");
  14.     z3 = zsub(z1,z2);
  15.     if ( fabs(z1.re-z2.re-z3.re) + fabs(z1.im-z2.im-z3.im) > 10*MACHEPS )
  16.     errmesg("zadd");
  17.     z3 = zmlt(z1,z2);
  18.     if ( fabs(z1.re*z2.re - z1.im*z2.im - z3.re) +
  19.      fabs(z1.im*z2.re + z1.re*z2.im - z3.im) > 10*MACHEPS )
  20.     errmesg("zmlt");
  21.     s1 = zabs(z1);
  22.     if ( fabs(s1*s1 - (z1.re*z1.re+z1.im*z1.im)) > 10*MACHEPS )
  23.     errmesg("zabs");
  24.     if ( zabs(zsub(z1,zmlt(z2,zdiv(z1,z2)))) > 10*MACHEPS ||
  25.      zabs(zsub(ONE,zdiv(z1,zmlt(z2,zdiv(z1,z2))))) > 10*MACHEPS )
  26.     errmesg("zdiv");
  27.  
  28.     z3 = zsqrt(z1);
  29.     if ( zabs(zsub(z1,zmlt(z3,z3))) > 10*MACHEPS )
  30.     errmesg("zsqrt");
  31.     if ( zabs(zsub(z1,zlog(zexp(z1)))) > 10*MACHEPS )
  32.     errmesg("zexp/zlog");
  33.     
  34.  
  35.     printf("# Check: MACHEPS = %g\n",MACHEPS);
  36.     /* allocate, initialise, copy and resize operations */
  37.     /* ZVEC */
  38.     notice("vector initialise, copy & resize");
  39.     x = zv_get(12);
  40.     y = zv_get(15);
  41.     z = zv_get(12);
  42.     zv_rand(x);
  43.     zv_rand(y);
  44.     z = zv_copy(x,z);
  45.     if ( zv_norm2(zv_sub(x,z,z)) >= MACHEPS )
  46.     errmesg("ZVEC copy");
  47.     zv_copy(x,y);
  48.     x = zv_resize(x,10);
  49.     y = zv_resize(y,10);
  50.     if ( zv_norm2(zv_sub(x,y,z)) >= MACHEPS )
  51.     errmesg("ZVEC copy/resize");
  52.     x = zv_resize(x,15);
  53.     y = zv_resize(y,15);
  54.     if ( zv_norm2(zv_sub(x,y,z)) >= MACHEPS )
  55.     errmesg("VZEC resize");
  56.  
  57.     /* ZMAT */
  58.     notice("matrix initialise, copy & resize");
  59.     A = zm_get(8,5);
  60.     B = zm_get(3,9);
  61.     C = zm_get(8,5);
  62.     zm_rand(A);
  63.     zm_rand(B);
  64.     C = zm_copy(A,C);
  65.     if ( zm_norm_inf(zm_sub(A,C,C)) >= MACHEPS )
  66.     errmesg("ZMAT copy");
  67.     zm_copy(A,B);
  68.     A = zm_resize(A,3,5);
  69.     B = zm_resize(B,3,5);
  70.     if ( zm_norm_inf(zm_sub(A,B,C)) >= MACHEPS )
  71.     errmesg("ZMAT copy/resize");
  72.     A = zm_resize(A,10,10);
  73.     B = zm_resize(B,10,10);
  74.     if ( zm_norm_inf(zm_sub(A,B,C)) >= MACHEPS )
  75.     errmesg("ZMAT resize");
  76.  
  77.     MEMCHK();
  78.  
  79.     /* PERM */
  80.     notice("permutation initialise, inverting & permuting vectors");
  81.     pi1 = px_get(15);
  82.     pi2 = px_get(12);
  83.     px_rand(pi1);
  84.     zv_rand(x);
  85.     px_zvec(pi1,x,z);
  86.     y = zv_resize(y,x->dim);
  87.     pxinv_zvec(pi1,z,y);
  88.     if ( zv_norm2(zv_sub(x,y,z)) >= MACHEPS )
  89.     errmesg("PERMute vector");
  90.  
  91.     /* testing catch() etc */
  92.     notice("error handling routines");
  93.     catch(E_NULL,
  94.       catchall(zv_add(ZVNULL,ZVNULL,ZVNULL);
  95.              errmesg("tracecatch() failure"),
  96.              printf("# tracecatch() caught error\n");
  97.              error(E_NULL,"main"));
  98.                  errmesg("catch() failure"),
  99.       printf("# catch() caught E_NULL error\n"));
  100.  
  101.     /* testing inner products and v_mltadd() etc */
  102.     notice("inner products and linear combinations");
  103.     u = zv_get(x->dim);
  104.     zv_rand(u);
  105.     zv_rand(x);
  106.     zv_resize(y,x->dim);
  107.     zv_rand(y);
  108.     zv_mltadd(y,x,zneg(zdiv(zin_prod(x,y),zin_prod(x,x))),z);
  109.     if ( zabs(zin_prod(x,z)) >= 5*MACHEPS*x->dim )
  110.     {
  111.     errmesg("zv_mltadd()/zin_prod()");
  112.     printf("# error norm = %g\n", zabs(zin_prod(x,z)));
  113.     }
  114.  
  115.     z1 = zneg(zdiv(zin_prod(x,y),zmake(zv_norm2(x)*zv_norm2(x),0.0)));
  116.     zv_mlt(z1,x,u);
  117.     zv_add(y,u,u);
  118.     if ( zv_norm2(zv_sub(u,z,u)) >= MACHEPS*x->dim )
  119.     {
  120.     errmesg("zv_mlt()/zv_norm2()");
  121.     printf("# error norm = %g\n", zv_norm2(u));
  122.     }
  123.  
  124. #ifdef ANSI_C
  125.     zv_linlist(u,x,z1,y,ONE,VNULL);
  126.     if ( zv_norm2(zv_sub(u,z,u)) >= MACHEPS*x->dim )
  127.     errmesg("zv_linlist()");
  128. #endif
  129. #ifdef VARARGS
  130.     zv_linlist(u,x,z1,y,ONE,VNULL);
  131.     if ( zv_norm2(zv_sub(u,z,u)) >= MACHEPS*x->dim )
  132.     errmesg("zv_linlist()");
  133. #endif
  134.  
  135.     MEMCHK();
  136.  
  137.     /* vector norms */
  138.     notice("vector norms");
  139.     x = zv_resize(x,12);
  140.     zv_rand(x);
  141.     for ( i = 0; i < x->dim; i++ )
  142.     if ( zabs(v_entry(x,i)) >= 0.7 )
  143.         v_set_val(x,i,ONE);
  144.         else
  145.         v_set_val(x,i,zneg(ONE));
  146.     s1 = zv_norm1(x);
  147.     s2 = zv_norm2(x);    
  148.     s3 = zv_norm_inf(x);
  149.     if ( fabs(s1 - x->dim) >= MACHEPS*x->dim ||
  150.      fabs(s2 - sqrt((double)(x->dim))) >= MACHEPS*x->dim ||
  151.      fabs(s3 - 1.0) >= MACHEPS )
  152.     errmesg("zv_norm1/2/_inf()");
  153.  
  154.     /* test matrix multiply etc */
  155.     notice("matrix multiply and invert");
  156.     A = zm_resize(A,10,10);
  157.     B = zm_resize(B,10,10);
  158.     zm_rand(A);
  159.     zm_inverse(A,B);
  160.     zm_mlt(A,B,C);
  161.     for ( i = 0; i < C->m; i++ )
  162.     m_set_val(C,i,i,zsub(m_entry(C,i,i),ONE));
  163.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  164.     errmesg("zm_inverse()/zm_mlt()");
  165.  
  166.     MEMCHK();
  167.  
  168.     /* ... and adjoints */
  169.     notice("adjoints and adjoint-multiplies");
  170.     zm_adjoint(A,A);    /* can do square matrices in situ */
  171.     zmam_mlt(A,B,C);
  172.     for ( i = 0; i < C->m; i++ )
  173.     m_set_val(C,i,i,zsub(m_entry(C,i,i),ONE));
  174.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  175.     errmesg("zm_adjoint()/zmam_mlt()");
  176.     zm_adjoint(A,A);
  177.     zm_adjoint(B,B);
  178.     zmma_mlt(A,B,C);
  179.     for ( i = 0; i < C->m; i++ )
  180.     m_set_val(C,i,i,zsub(m_entry(C,i,i),ONE));
  181.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  182.     errmesg("zm_adjoint()/zmma_mlt()");
  183.     zsm_mlt(zmake(3.71,2.753),B,B);
  184.     zmma_mlt(A,B,C);
  185.     for ( i = 0; i < C->m; i++ )
  186.     m_set_val(C,i,i,zsub(m_entry(C,i,i),zmake(3.71,-2.753)));
  187.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  188.     errmesg("szm_mlt()/zmma_mlt()");
  189.     zm_adjoint(B,B);
  190.     zsm_mlt(zdiv(ONE,zmake(3.71,-2.753)),B,B);
  191.  
  192.     MEMCHK();
  193.  
  194.     /* ... and matrix-vector multiplies */
  195.     notice("matrix-vector multiplies");
  196.     x = zv_resize(x,A->n);
  197.     y = zv_resize(y,A->m);
  198.     z = zv_resize(z,A->m);
  199.     u = zv_resize(u,A->n);
  200.     zv_rand(x);
  201.     zv_rand(y);
  202.     zmv_mlt(A,x,z);
  203.     z1 = zin_prod(y,z);
  204.     zvm_mlt(A,y,u);
  205.     z2 = zin_prod(u,x);
  206.     if ( zabs(zsub(z1,z2)) >= (MACHEPS*x->dim)*x->dim )
  207.     {
  208.     errmesg("zmv_mlt()/zvm_mlt()");
  209.     printf("# difference between inner products is %g\n",
  210.            zabs(zsub(z1,z2)));
  211.     }
  212.     zmv_mlt(B,z,u);
  213.     if ( zv_norm2(zv_sub(u,x,u)) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  214.     errmesg("zmv_mlt()/zvm_mlt()");
  215.  
  216.     MEMCHK();
  217.  
  218.     /* get/set row/col */
  219.     notice("getting and setting rows and cols");
  220.     x = zv_resize(x,A->n);
  221.     y = zv_resize(y,B->m);
  222.     x = zget_row(A,3,x);
  223.     y = zget_col(B,3,y);
  224.     if ( zabs(zsub(_zin_prod(x,y,0,Z_NOCONJ),ONE)) >=
  225.     MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  226.     errmesg("zget_row()/zget_col()");
  227.     zv_mlt(zmake(-1.0,0.0),x,x);
  228.     zv_mlt(zmake(-1.0,0.0),y,y);
  229.     zset_row(A,3,x);
  230.     zset_col(B,3,y);
  231.     zm_mlt(A,B,C);
  232.     for ( i = 0; i < C->m; i++ )
  233.     m_set_val(C,i,i,zsub(m_entry(C,i,i),ONE));
  234.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  235.     errmesg("zset_row()/zset_col()");
  236.  
  237.     MEMCHK();
  238.  
  239.     /* matrix norms */
  240.     notice("matrix norms");
  241.     A = zm_resize(A,11,15);
  242.     zm_rand(A);
  243.     s1 = zm_norm_inf(A);
  244.     B = zm_adjoint(A,B);
  245.     s2 = zm_norm1(B);
  246.     if ( fabs(s1 - s2) >= MACHEPS*A->m )
  247.     errmesg("zm_norm1()/zm_norm_inf()");
  248.     C = zmam_mlt(A,A,C);
  249.     z1.re = z1.im = 0.0;
  250.     for ( i = 0; i < C->m && i < C->n; i++ )
  251.     z1 = zadd(z1,m_entry(C,i,i));
  252.     if ( fabs(sqrt(z1.re) - zm_norm_frob(A)) >= MACHEPS*A->m*A->n )
  253.     errmesg("zm_norm_frob");
  254.  
  255.     MEMCHK();
  256.     
  257.     /* permuting rows and columns */
  258.     /******************************
  259.     notice("permuting rows & cols");
  260.     A = zm_resize(A,11,15);
  261.     B = zm_resize(B,11,15);
  262.     pi1 = px_resize(pi1,A->m);
  263.     px_rand(pi1);
  264.     x = zv_resize(x,A->n);
  265.     y = zmv_mlt(A,x,y);
  266.     px_rows(pi1,A,B);
  267.     px_zvec(pi1,y,z);
  268.     zmv_mlt(B,x,u);
  269.     if ( zv_norm2(zv_sub(z,u,u)) >= MACHEPS*A->m )
  270.     errmesg("px_rows()");
  271.     pi1 = px_resize(pi1,A->n);
  272.     px_rand(pi1);
  273.     px_cols(pi1,A,B);
  274.     pxinv_zvec(pi1,x,z);
  275.     zmv_mlt(B,z,u);
  276.     if ( zv_norm2(zv_sub(y,u,u)) >= MACHEPS*A->n )
  277.     errmesg("px_cols()");
  278.     ******************************/
  279.  
  280.     MEMCHK();
  281.  
  282.     /* MATLAB save/load */
  283.     notice("MATLAB save/load");
  284.     A = zm_resize(A,12,11);
  285.     if ( (fp=fopen(SAVE_FILE,"w")) == (FILE *)NULL )
  286.     printf("Cannot perform MATLAB save/load test\n");
  287.     else
  288.     {
  289.     zm_rand(A);
  290.     zm_save(fp, A, name);
  291.     fclose(fp);
  292.     if ( (fp=fopen(SAVE_FILE,"r")) == (FILE *)NULL )
  293.         printf("Cannot open save file \"%s\"\n",SAVE_FILE);
  294.     else
  295.     {
  296.         ZM_FREE(B);
  297.         B = zm_load(fp,&cp);
  298.         if ( strcmp(name,cp) || zm_norm1(zm_sub(A,B,C)) >=
  299.          MACHEPS*A->m )
  300.         {
  301.         errmesg("zm_load()/zm_save()");
  302.         printf("# orig. name = %s, restored name = %s\n", name, cp);
  303.         printf("# orig. A =\n");    zm_output(A);
  304.         printf("# restored A =\n");    zm_output(B);
  305.         }
  306.     }
  307.     }
  308.  
  309.     MEMCHK();
  310.  
  311.     /* Now, onto matrix factorisations */
  312.     A = zm_resize(A,10,10);
  313.     B = zm_resize(B,A->m,A->n);
  314.     zm_copy(A,B);
  315.     x = zv_resize(x,A->n);
  316.     y = zv_resize(y,A->m);
  317.     z = zv_resize(z,A->n);
  318.     u = zv_resize(u,A->m);
  319.     zv_rand(x);
  320.     zmv_mlt(B,x,y);
  321.     z = zv_copy(x,z);
  322.  
  323.     notice("LU factor/solve");
  324.     pivot = px_get(A->m);
  325.     zLUfactor(A,pivot);
  326.     tracecatch(zLUsolve(A,pivot,y,x),"main");
  327.     tracecatch(cond_est = zLUcondest(A,pivot),"main");
  328.     printf("# cond(A) approx= %g\n", cond_est);
  329.     if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est)
  330.     {
  331.     errmesg("zLUfactor()/zLUsolve()");
  332.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  333.            zv_norm2(zv_sub(x,z,u)), MACHEPS);
  334.     }
  335.  
  336.  
  337.     zv_copy(y,x);
  338.     tracecatch(zLUsolve(A,pivot,x,x),"main");
  339.     tracecatch(cond_est = zLUcondest(A,pivot),"main");
  340.     if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est)
  341.     {
  342.     errmesg("zLUfactor()/zLUsolve()");
  343.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  344.            zv_norm2(zv_sub(x,z,u)), MACHEPS);
  345.     }
  346.  
  347.     zvm_mlt(B,z,y);
  348.     zv_copy(y,x);
  349.     tracecatch(zLUAsolve(A,pivot,x,x),"main");
  350.     if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est)
  351.     {
  352.     errmesg("zLUfactor()/zLUAsolve()");
  353.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  354.            zv_norm2(zv_sub(x,z,u)), MACHEPS);
  355.     }
  356.  
  357.     MEMCHK();
  358.  
  359.     /* QR factorisation */
  360.     zm_copy(B,A);
  361.     zmv_mlt(B,z,y);
  362.     notice("QR factor/solve:");
  363.     diag = zv_get(A->m);
  364.     zQRfactor(A,diag);
  365.     zQRsolve(A,diag,y,x);
  366.     if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est )
  367.     {
  368.     errmesg("zQRfactor()/zQRsolve()");
  369.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  370.            zv_norm2(zv_sub(x,z,u)), MACHEPS);
  371.     }
  372.     printf("# QR cond(A) approx= %g\n", zQRcondest(A));
  373.     Q = zm_get(A->m,A->m);
  374.     zmakeQ(A,diag,Q);
  375.     zmakeR(A,A);
  376.     zm_mlt(Q,A,C);
  377.     zm_sub(B,C,C);
  378.     if ( zm_norm1(C) >= MACHEPS*zm_norm1(Q)*zm_norm1(B) )
  379.     {
  380.     errmesg("zQRfactor()/zmakeQ()/zmakeR()");
  381.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  382.            zm_norm1(C), MACHEPS);
  383.     }
  384.  
  385.     MEMCHK();
  386.  
  387.     /* now try with a non-square matrix */
  388.     A = zm_resize(A,15,7);
  389.     zm_rand(A);
  390.     B = zm_copy(A,B);
  391.     diag = zv_resize(diag,A->n);
  392.     x = zv_resize(x,A->n);
  393.     y = zv_resize(y,A->m);
  394.     zv_rand(y);
  395.     zQRfactor(A,diag);
  396.     x = zQRsolve(A,diag,y,x);
  397.     /* z is the residual vector */
  398.     zmv_mlt(B,x,z);    zv_sub(z,y,z);
  399.     /* check B*.z = 0 */
  400.     zvm_mlt(B,z,u);
  401.     if ( zv_norm2(u) >= 100*MACHEPS*zm_norm1(B)*zv_norm2(y) )
  402.     {
  403.     errmesg("zQRfactor()/zQRsolve()");
  404.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  405.            zv_norm2(u), MACHEPS);
  406.     }
  407.     Q = zm_resize(Q,A->m,A->m);
  408.     zmakeQ(A,diag,Q);
  409.     zmakeR(A,A);
  410.     zm_mlt(Q,A,C);
  411.     zm_sub(B,C,C);
  412.     if ( zm_norm1(C) >= MACHEPS*zm_norm1(Q)*zm_norm1(B) )
  413.     {
  414.     errmesg("zQRfactor()/zmakeQ()/zmakeR()");
  415.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  416.            zm_norm1(C), MACHEPS);
  417.     }
  418.     D = zm_get(A->m,Q->m);
  419.     zmam_mlt(Q,Q,D);
  420.     for ( i = 0; i < D->m; i++ )
  421.     m_set_val(D,i,i,zsub(m_entry(D,i,i),ONE));
  422.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q) )
  423.     {
  424.     errmesg("QRfactor()/makeQ()/makeR()");
  425.     printf("# QR orthogonality error = %g [cf MACHEPS = %g]\n",
  426.            zm_norm1(D), MACHEPS);
  427.     }
  428.  
  429.     MEMCHK();
  430.  
  431.     /* QRCP factorisation */
  432.     zm_copy(B,A);
  433.     notice("QR factor/solve with column pivoting");
  434.     pivot = px_resize(pivot,A->n);
  435.     zQRCPfactor(A,diag,pivot);
  436.     z = zv_resize(z,A->n);
  437.     zQRCPsolve(A,diag,pivot,y,z);
  438.     /* pxinv_zvec(pivot,z,x); */
  439.     /* now compute residual (z) vector */
  440.     zmv_mlt(B,x,z);    zv_sub(z,y,z);
  441.     /* check B^T.z = 0 */
  442.     zvm_mlt(B,z,u);
  443.     if ( zv_norm2(u) >= MACHEPS*zm_norm1(B)*zv_norm2(y) )
  444.     {
  445.     errmesg("QRCPfactor()/QRsolve()");
  446.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  447.            zv_norm2(u), MACHEPS);
  448.     }
  449.  
  450.     Q = zm_resize(Q,A->m,A->m);
  451.     zmakeQ(A,diag,Q);
  452.     zmakeR(A,A);
  453.     zm_mlt(Q,A,C);
  454.     ZM_FREE(D);
  455.     D = zm_get(B->m,B->n);
  456.     /******************************
  457.     px_cols(pivot,C,D);
  458.     zm_sub(B,D,D);
  459.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm1(B) )
  460.     {
  461.     errmesg("QRCPfactor()/makeQ()/makeR()");
  462.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  463.            zm_norm1(D), MACHEPS);
  464.     }
  465.     ******************************/
  466.  
  467.     /* Now check eigenvalue/SVD routines */
  468.     notice("complex Schur routines");
  469.     A = zm_resize(A,11,11);
  470.     B = zm_resize(B,A->m,A->n);
  471.     C = zm_resize(C,A->m,A->n);
  472.     D = zm_resize(D,A->m,A->n);
  473.     Q = zm_resize(Q,A->m,A->n);
  474.  
  475.     MEMCHK();
  476.  
  477.     /* now test complex Schur decomposition */
  478.     /* zm_copy(A,B); */
  479.     ZM_FREE(A);
  480.     A = zm_get(11,11);
  481.     zm_rand(A);
  482.     B = zm_copy(A,B);
  483.     MEMCHK();
  484.  
  485.     B = zschur(B,Q);
  486.     MEMCHK();
  487.  
  488.     zm_mlt(Q,B,C);
  489.     zmma_mlt(C,Q,D);
  490.     MEMCHK();
  491.     zm_sub(A,D,D);
  492.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q)*zm_norm1(B)*5 )
  493.     {
  494.     errmesg("zschur()");
  495.     printf("# Schur reconstruction error = %g [cf MACHEPS = %g]\n",
  496.            zm_norm1(D), MACHEPS);
  497.     }
  498.  
  499.     /* orthogonality check */
  500.     zmma_mlt(Q,Q,D);
  501.     for ( i = 0; i < D->m; i++ )
  502.     m_set_val(D,i,i,zsub(m_entry(D,i,i),ONE));
  503.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q)*10 )
  504.     {
  505.     errmesg("zschur()");
  506.     printf("# Schur orthogonality error = %g [cf MACHEPS = %g]\n",
  507.            zm_norm1(D), MACHEPS);
  508.     }
  509.  
  510.     MEMCHK();
  511.  
  512.     /* now test SVD */
  513.     /******************************
  514.     A = zm_resize(A,11,7);
  515.     zm_rand(A);
  516.     U = zm_get(A->n,A->n);
  517.     Q = zm_resize(Q,A->m,A->m);
  518.     u = zv_resize(u,max(A->m,A->n));
  519.     svd(A,Q,U,u);
  520.     ******************************/
  521.     /* check reconstruction of A */
  522.     /******************************
  523.     D = zm_resize(D,A->m,A->n);
  524.     C = zm_resize(C,A->m,A->n);
  525.     zm_zero(D);
  526.     for ( i = 0; i < min(A->m,A->n); i++ )
  527.     zm_set_val(D,i,i,v_entry(u,i));
  528.     zmam_mlt(Q,D,C);
  529.     zm_mlt(C,U,D);
  530.     zm_sub(A,D,D);
  531.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(U)*zm_norm_inf(Q)*zm_norm1(A) )
  532.     {
  533.     errmesg("svd()");
  534.     printf("# SVD reconstruction error = %g [cf MACHEPS = %g]\n",
  535.            zm_norm1(D), MACHEPS);
  536.     }
  537.     ******************************/
  538.     /* check orthogonality of Q and U */
  539.     /******************************
  540.     D = zm_resize(D,Q->n,Q->n);
  541.     zmam_mlt(Q,Q,D);
  542.     for ( i = 0; i < D->m; i++ )
  543.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  544.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q)*5 )
  545.     {
  546.     errmesg("svd()");
  547.     printf("# SVD orthognality error (Q) = %g [cf MACHEPS = %g\n",
  548.            zm_norm1(D), MACHEPS);
  549.     }
  550.     D = zm_resize(D,U->n,U->n);
  551.     zmam_mlt(U,U,D);
  552.     for ( i = 0; i < D->m; i++ )
  553.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  554.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(U)*zm_norm_inf(U)*5 )
  555.     {
  556.     errmesg("svd()");
  557.     printf("# SVD orthognality error (U) = %g [cf MACHEPS = %g\n",
  558.            zm_norm1(D), MACHEPS);
  559.     }
  560.     for ( i = 0; i < u->dim; i++ )
  561.     if ( v_entry(u,i) < 0 || (i < u->dim-1 &&
  562.                   v_entry(u,i+1) > v_entry(u,i)) )
  563.         break;
  564.     if ( i < u->dim )
  565.     {
  566.     errmesg("svd()");
  567.     printf("# SVD sorting error\n");
  568.     }
  569.     ******************************/
  570.  
  571.     ZV_FREE(x);    ZV_FREE(y);    ZV_FREE(z);
  572.     ZV_FREE(u);    ZV_FREE(diag);
  573.     PX_FREE(pi1);    PX_FREE(pi2);    PX_FREE(pivot);
  574.     ZM_FREE(A);    ZM_FREE(B);    ZM_FREE(C);
  575.     ZM_FREE(D);    ZM_FREE(Q);
  576.  
  577.     mem_stat_free(1);
  578.  
  579.     MEMCHK();
  580.     printf("# Finished torture test for complex numbers/vectors/matrices\n");
  581.     mem_info();
  582. }
  583.  
  584. FileDataŵzvecopÙ+Eÿÿÿ(”3¼1
  585. /**************************************************************************
  586. **
  587. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  588. **
  589. **                 Meschach Library
  590. ** 
  591. ** This Meschach Library is provided "as is" without any express 
  592. ** or implied warranty of any kind with respect to this software. 
  593. ** In particular the authors shall not be liable for any direct, 
  594. ** indirect, special, incidental or consequential damages arising 
  595. ** in any way from use of the software.
  596. ** 
  597. ** Everyone is granted permission to copy, modify and redistribute this
  598. ** Meschach Library, provided:
  599. **  1.  All copies contain this copyright notice.
  600. **  2.  All modified copies shall carry a notice stating who
  601. **      made the last modification and the date of such modification.
  602. **  3.  No charge is made for this software or works derived from it.  
  603. **      This clause shall not be construed as constraining other software
  604. **