home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / zqrfctr < prev    next >
Encoding:
Text File  |  1994-08-04  |  13.5 KB  |  507 lines
  1.  0.0 )
  2.         sum.re = 1.0;
  3.     else
  4.     {
  5.         sum.re += sum.re / norm1;
  6.         sum.im += sum.im / norm1;
  7.     }
  8.     /* y->ve[i] = sum / QR->me[i][i]; */
  9.     y->ve[i] = zdiv(sum,QR->me[i][i]);
  10.     }
  11.     zUAmlt(QR,y,y);
  12.  
  13.     /* no* indirect, special, incidental or consequential damages arising 
  14. ** in any way from use of the software.
  15. ** 
  16. ** Everyone is granted permission to copy, modify and redistribute this
  17. ** Meschach Library, provided:
  18. **  1.  All copies contain this copyright notice.
  19. **  2.  All modified copies shall carry a notice stating who
  20. **      made the last modification and the date of such modification.
  21. **  3.  No charge is made for this software or works derived from it.  
  22. **      This clause shall not be construed as constraining other software
  23. **      distributed on the same medium as this software, nor is a
  24. **      distribution fee considered a charge.
  25. **
  26. ***************************************************************************/
  27.  
  28.  
  29. /*
  30.     This file contains routines for import/exporting complex data
  31.     to/from MATLAB. The main routines are:
  32.             ZMAT *zm_save(FILE *fp,ZMAT *A,char *name)
  33.             ZVEC *zv_save(FILE *fp,ZVEC *x,char *name)
  34.             complex z_save(FILE *fp,complex z,char *name)
  35.             ZMAT *zm_load(FILE *fp,char **name)
  36. */
  37.  
  38. #include        <stdio.h>
  39. #include        "zmatrix.h"
  40. #include    "matlab.h"
  41.  
  42. static char rcsid[] = "$Id: zmatlab.c,v 1.1 1994/01/13 04:24:57 des Exp $";
  43.  
  44. /* zm_save -- save matrix in ".mat" file for MATLAB
  45.    -- returns matrix to be saved */
  46. ZMAT    *zm_save(fp,A,name)
  47. FILE    *fp;
  48. ZMAT    *A;
  49. char    *name;
  50. {
  51.     int     i, j;
  52.     matlab  mat;
  53.     
  54.     if ( ! A )
  55.     error(E_NULL,"zm_save");
  56.     
  57.     mat.type = 1000*MACH_ID + 100*ORDER + 10*PRECISION + 0;
  58.     mat.m = A->m;
  59.     mat.n = A->n;
  60.     mat.imag = TRUE;
  61.     mat.namlen = (name == (char *)NULL) ? 1 : strlen(name)+1;
  62.     
  63.     /* write header */
  64.     fwrite(&mat,sizeof(matlab),1,fp);
  65.     /* write name */
  66.     if ( name == (char *)NULL )
  67.     fwrite("",sizeof(char),1,fp);
  68.     else
  69.     fwrite(name,sizeof(char),(int)(mat.namlen),fp);
  70.     /* write actual data */
  71.     for ( i = 0; i < A->m; i++ )
  72.     for ( j = 0; j < A->n; j++ )
  73.         fwrite(&(A->me[i][j].re),sizeof(Real),1,fp);
  74.     for ( i = 0; i < A->m; i++ )
  75.     for ( j = 0; j < A->n; j++ )
  76.         fwrite(&(A->me[i][j].im),sizeof(Real),1,fp);
  77.     
  78.     return A;
  79. }
  80.  
  81.  
  82. /* zv_save -- save vector in ".mat" file for MATLAB
  83.    -- saves it as a row vector
  84.    -- returns vector to be saved */
  85. ZVEC    *zv_save(fp,x,name)
  86. FILE    *fp;
  87. ZVEC    *x;
  88. char    *name;
  89. {
  90.     int    i;
  91.     matlab  mat;
  92.     
  93.     if ( ! x )
  94.     error(E_NULL,"zv_save");
  95.     
  96.     mat.type = 1000*MACH_ID + 100*ORDER + 10*PRECISION + 0;
  97.     mat.m = x->dim;
  98.     mat.n = 1;
  99.     mat.imag = TRUE;
  100.     mat.namlen = (name == (char *)NULL) ? 1 : strlen(name)+1;
  101.     
  102.     /* write header */
  103.     fwrite(&mat,sizeof(matlab),1,fp);
  104.     /* write name */
  105.     if ( name == (char *)NULL )
  106.     fwrite("",sizeof(char),1,fp);
  107.     else
  108.     fwrite(name,sizeof(char),(int)(mat.namlen),fp);
  109.     /* write actual data */
  110.     for ( i = 0; i < x->dim; i++ )
  111.     fwrite(&(x->ve[i].re),sizeof(Real),1,fp);
  112.     for ( i = 0; i < x->dim; i++ )
  113.     fwrite(&(x->ve[i].im),sizeof(Real),1,fp);
  114.     
  115.     return x;
  116. }
  117.  
  118. /* z_save -- saves complex number in ".mat" file for MATLAB
  119.     -- returns complex number to be saved */
  120. complex    z_save(fp,z,name)
  121. FILE    *fp;
  122. complex    z;
  123. char    *name;
  124. {
  125.     matlab  mat;
  126.     
  127.     mat.type = 1000*MACH_ID + 100*ORDER + 10*PRECISION + 0;
  128.     mat.m = 1;
  129.     mat.n = 1;
  130.     mat.imag = TRUE;
  131.     mat.namlen = (name == (char *)NULL) ? 1 : strlen(name)+1;
  132.     
  133.     /* write header */
  134.     fwrite(&mat,sizeof(matlab),1,fp);
  135.     /* write name */
  136.     if ( name == (char *)NULL )
  137.     fwrite("",sizeof(char),1,fp);
  138.     else
  139.     fwrite(name,sizeof(char),(int)(mat.namlen),fp);
  140.     /* write actual data */
  141.     fwrite(&z,sizeof(complex),1,fp);
  142.     
  143.     return z;
  144. }
  145.  
  146.  
  147.  
  148. /* zm_load -- loads in a ".mat" file variable as produced by MATLAB
  149.    -- matrix returned; imaginary parts ignored */
  150. ZMAT    *zm_load(fp,name)
  151. FILE    *fp;
  152. char    **name;
  153. {
  154.     ZMAT     *A;
  155.     int     i;
  156.     int     m_flag, o_flag, p_flag, t_flag;
  157.     float   f_temp;
  158.     double  d_temp;
  159.     matlab  mat;
  160.     
  161.     if ( fread(&mat,sizeof(matlab),1,fp) != 1 )
  162.     error(E_FORMAT,"zm_load");
  163.     if ( mat.type >= 10000 )    /* don't load a sparse matrix! */
  164.     error(E_FORMAT,"zm_load");
  165.     m_flag = (mat.type/1000) % 10;
  166.     o_flag = (mat.type/100) % 10;
  167.     p_flag = (mat.type/10) % 10;
  168.     t_flag = (mat.type) % 10;
  169.     if ( m_flag != MACH_ID )
  170.     error(E_FORMAT,"zm_load");
  171.     if ( t_flag != 0 )
  172.     error(E_FORMAT,"zm_load");
  173.     if ( p_flag != DOUBLE_PREC && p_flag != SINGLE_PREC )
  174.     error(E_FORMAT,"zm_load");
  175.     *name = (char *)malloc((unsigned)(mat.namlen)+1);
  176.     if ( fread(*name,sizeof(char),(unsigned)(mat.namlen),fp) == 0 )
  177.     error(E_FORMAT,"zm_load");
  178.     A = zm_get((unsigned)(mat.m),(unsigned)(mat.n));
  179.     for ( i = 0; i < A->m*A->n; i++ )
  180.     {
  181.     if ( p_flag == DOUBLE_PREC )
  182.         fread(&d_temp,sizeof(double),1,fp);
  183.     else
  184.     {
  185.         fread(&f_temp,sizeof(float),1,fp);
  186.         d_temp = f_temp;
  187.     }
  188.     if ( o_flag == ROW_ORDER )
  189.         A->me[i / A->n][i % A->n].re = d_temp;
  190.     else if ( o_flag == COL_ORDER )
  191.         A->me[i % A->m][i / A->m].re = d_temp;
  192.     else
  193.         error(E_FORMAT,"zm_load");
  194.     }
  195.     
  196.     if ( mat.imag )         /* skip imaginary part */
  197.     for ( i = 0; i < A->m*A->n; i++ )
  198.     {
  199.         if ( p_flag == DOUBLE_PREC )
  200.         fread(&d_temp,sizeof(double),1,fp);
  201.         else
  202.         {
  203.         fread(&f_temp,sizeof(float),1,fp);
  204.         d_temp = f_temp;
  205.         }
  206.         if ( o_flag == ROW_ORDER )
  207.         A->me[i / A->n][i % A->n].im = d_temp;
  208.         else if ( o_flag == COL_ORDER )
  209.         A->me[i % A->m][i / A->m].im = d_temp;
  210.         else
  211.         error(E_FORMAT,"zm_load");
  212.     }
  213.     
  214.     return A;
  215. }
  216.  
  217. FileDataŵzmatop”;Eýÿÿ õ 
  218. /**************************************************************************
  219. **
  220. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  221. **
  222. **                 Meschach Library
  223. ** 
  224. ** This Meschach Library is provided "as is" without any express 
  225. ** or implied warranty of any kind with respect to this software. 
  226. ** In particular the authors shall not be liable for any direct, 
  227. ** indirect, special, incidental or consequential damages arising 
  228. ** in any way from use of the software.
  229. ** 
  230. ** Everyone is granted permission to copy, modify and redistribute this
  231. ** Meschach Library, provided:
  232. **  1.  All copies contain this copyright notice.
  233. **  2.  All modified copies shall carry a notice stating who
  234. **      made the last modification and the date of such modification.
  235. **  3.  No charge is made for this software or works derived from it.  
  236. **      This clause shall not be construed as constraining other software
  237. **      distributed on the same medium as this software, nor is a
  238. **      distribution fee considered a charge.
  239. **
  240. ***************************************************************************/
  241.  
  242.  
  243.  
  244. #include    <stdio.h>
  245. #include    "zmatrix.h"
  246.  
  247. static    char    rcsid[] = "$Id: zmatop.c,v 1.1 1994/01/13 04:24:46 des Exp $";
  248.  
  249.  
  250. #define    is_zero(z)    ((z).re == 0.0 && (z).im == 0.0)
  251.  
  252. /* zm_add -- matrix addition -- may be in-situ */
  253. ZMAT    *zm_add(mat1,mat2,out)
  254. ZMAT    *mat1,*mat2,*out;
  255. {
  256.     u_int    m,n,i;
  257.     
  258.     if ( mat1==ZMNULL || mat2==ZMNULL )
  259.     error(E_NULL,"zm_add");
  260.     if ( mat1->m != mat2->m || mat1->n != mat2->n )
  261.     error(E_SIZES,"zm_add");
  262.     if ( out==ZMNULL || out->m != mat1->m || out->n != mat1->n )
  263.     out = zm_resize(out,mat1->m,mat1->n);
  264.     m = mat1->m;    n = mat1->n;
  265.     for ( i=0; i<m; i++ )
  266.     {
  267.     __zadd__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  268.     /**************************************************
  269.       for ( j=0; j<n; j++ )
  270.       out->me[i][j] = mat1->me[i][j]+mat2->me[i][j];
  271.       **************************************************/
  272.     }
  273.     
  274.     return (out);
  275. }
  276.  
  277. /* zm_sub -- matrix subtraction -- may be in-situ */
  278. ZMAT    *zm_sub(mat1,mat2,out)
  279. ZMAT    *mat1,*mat2,*out;
  280. {
  281.     u_int    m,n,i;
  282.     
  283.     if ( mat1==ZMNULL || mat2==ZMNULL )
  284.     error(E_NULL,"zm_sub");
  285.     if ( mat1->m != mat2->m || mat1->n != mat2->n )
  286.     error(E_SIZES,"zm_sub");
  287.     if ( out==ZMNULL || out->m != mat1->m || out->n != mat1->n )
  288.     out = zm_resize(out,mat1->m,mat1->n);
  289.     m = mat1->m;    n = mat1->n;
  290.     for ( i=0; i<m; i++ )
  291.     {
  292.     __zsub__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  293.     /**************************************************
  294.       for ( j=0; j<n; j++ )
  295.       out->me[i][j] = mat1->me[i][j]-mat2->me[i][j];
  296.     **************************************************/
  297.     }
  298.     
  299.     return (out);
  300. }
  301.  
  302. /*
  303.   Note: In the following routines, "adjoint" means complex conjugate
  304.   transpose:
  305.   A* = conjugate(A^T)
  306.   */
  307.  
  308. /* zm_mlt -- matrix-matrix multiplication */
  309. ZMAT    *zm_mlt(A,B,OUT)
  310. ZMAT    *A,*B,*OUT;
  311. {
  312.     u_int    i, /* j, */ k, m, n, p;
  313.     complex    **A_v, **B_v /*, *B_row, *OUT_row, sum, tmp */;
  314.     
  315.     if ( A==ZMNULL || B==ZMNULL )
  316.     error(E_NULL,"zm_mlt");
  317.     if ( A->n != B->m )
  318.     error(E_SIZES,"zm_mlt");
  319.     if ( A == OUT || B == OUT )
  320.     error(E_INSITU,"zm_mlt");
  321.     m = A->m;    n = A->n;    p = B->n;
  322.     A_v = A->me;        B_v = B->me;
  323.     
  324.     if ( OUT==ZMNULL || OUT->m != A->m || OUT->n != B->n )
  325.     OUT = zm_resize(OUT,A->m,B->n);
  326.     
  327.     /****************************************************************
  328.       for ( i=0; i<m; i++ )
  329.       for  ( j=0; j<p; j++ )
  330.       {
  331.       sum = 0.0;
  332.       for ( k=0; k<n; k++ )
  333.       sum += A_v[i][k]*B_v[k][j];
  334.       OUT->me[i][j] = sum;
  335.       }
  336.     ****************************************************************/
  337.     zm_zero(OUT);
  338.     for ( i=0; i<m; i++ )
  339.     for ( k=0; k<n; k++ )
  340.     {
  341.         if ( ! is_zero(A_v[i][k]) )
  342.         __zmltadd__(OUT->me[i],B_v[k],A_v[i][k],(int)p,Z_NOCONJ);
  343.         /**************************************************
  344.           B_row = B_v[k];    OUT_row = OUT->me[i];
  345.           for ( j=0; j<p; j++ )
  346.           (*OUT_row++) += tmp*(*B_row++);
  347.         **************************************************/
  348.     }
  349.     
  350.     return OUT;
  351. }
  352.  
  353. /* zmma_mlt -- matrix-matrix adjoint multiplication
  354.    -- A.B* is returned, and stored in OUT */
  355. ZMAT    *zmma_mlt(A,B,OUT)
  356. ZMAT    *A, *B, *OUT;
  357. {
  358.     int    i, j, limit;
  359.     /* complex    *A_row, *B_row, sum; */
  360.     
  361.     if ( ! A || ! B )
  362.     error(E_NULL,"zmma_mlt");
  363.     if ( A == OUT || B == OUT )
  364.     error(E_INSITU,"zmma_mlt");
  365.     if ( A->n != B->n )
  366.     error(E_SIZES,"zmma_mlt");
  367.     if ( ! OUT || OUT->m != A->m || OUT->n != B->m )
  368.     OUT = zm_resize(OUT,A->m,B->m);
  369.     
  370.     limit = A->n;
  371.     for ( i = 0; i < A->m; i++ )
  372.     for ( j = 0; j < B->m; j++ )
  373.     {
  374.         OUT->me[i][j] = __zip__(B->me[j],A->me[i],(int)limit,Z_CONJ);
  375.         /**************************************************
  376.           sum = 0.0;
  377.           A_row = A->me[i];
  378.           B_row = B->me[j];
  379.           for ( k = 0; k < limit; k++ )
  380.           sum += (*A_row++)*(*B_row++);
  381.           OUT->me[i][j] = sum;
  382.           **************************************************/
  383.     }
  384.     
  385.     return OUT;
  386. }
  387.  
  388. /* zmam_mlt -- matrix adjoint-matrix multiplication
  389.    -- A*.B is returned, result stored in OUT */
  390. ZMAT    *zmam_mlt(A,B,OUT)
  391. ZMAT    *A, *B, *OUT;
  392. {
  393.     int    i, k, limit;
  394.     /* complex    *B_row, *OUT_row, multiplier; */
  395.     complex    tmp;
  396.     
  397.     if ( ! A || ! B )
  398.     error(E_NULL,"zmam_mlt");
  399.     if ( A == OUT || B == OUT )
  400.     error(E_INSITU,"zmam_mlt");
  401.     if ( A->m != B->m )
  402.     error(E_SIZES,"zmam_mlt");
  403.     if ( ! OUT || OUT->m != A->n || OUT->n != B->n )
  404.     OUT = zm_resize(OUT,A->n,B->n);
  405.     
  406.     limit = B->n;
  407.     zm_zero(OUT);
  408.     for ( k = 0; k < A->m; k++ )
  409.     for ( i = 0; i < A->n; i++ )
  410.     {
  411.         tmp.re =   A->me[k][i].re;
  412.         tmp.im = - A->me[k][i].im;
  413.         if ( ! is_zero(tmp) )
  414.         __zmltadd__(OUT->me[i],B->me[k],tmp,(int)limit,Z_NOCONJ);
  415.     }
  416.     
  417.     return OUT;
  418. }
  419.  
  420. /* zmv_mlt -- matrix-vector multiplication 
  421.    -- Note: b is treated as a column vector */
  422. ZVEC    *zmv_mlt(A,b,out)
  423. ZMAT    *A;
  424. ZVEC    *b,*out;
  425. {
  426.     u_int    i, m, n;
  427.     complex    **A_v, *b_v /*, *A_row */;
  428.     /* register complex    sum; */
  429.     
  430.     if ( A==ZMNULL || b==ZVNULL )
  431.     error(E_NULL,"zmv_mlt");
  432.     if ( A->n != b->dim )
  433.     error(E_SIZES,"zmv_mlt");
  434.     if ( b == out )
  435.     error(E_INSITU,"zmv_mlt");
  436.     if ( out == ZVNULL || out->dim != A->m )
  437.     out = zv_resize(out,A->m);
  438.     
  439.     m = A->m;        n = A->n;
  440.     A_v = A->me;        b_v = b->ve;
  441.     for ( i=0; i<m; i++ )
  442.     {
  443.     /* for ( j=0; j<n; j++ )
  444.        sum += A_v[i][j]*b_v[j]; */
  445.     out->ve[i] = __zip__(A_v[i],b_v,(int)n,Z_NOCONJ);
  446.     /**************************************************
  447.       A_row = A_v[i];        b_v = b->ve;
  448.       for ( j=0; j<n; j++ )
  449.       sum += (*A_row++)*(*b_v++);
  450.       out->ve[i] = sum;
  451.     **************************************************/
  452.     }
  453.     
  454.     return out;
  455. }
  456.  
  457. /* zsm_mlt -- scalar-matrix multiply -- may be in-situ */
  458. ZMAT    *zsm_mlt(scalar,matrix,out)
  459. complex    scalar;
  460. ZMAT    *matrix,*out;
  461. {
  462.     u_int    m,n,i;
  463.     
  464.     if ( matrix==ZMNULL )
  465.     error(E_NULL,"zsm_mlt");
  466.     if ( out==ZMNULL || out->m != matrix->m || out->n != matrix->n )
  467.     out = zm_resize(out,matrix->m,matrix->n);
  468.     m = matrix->m;    n = matrix->n;
  469.     for ( i=0; i<m; i++ )
  470.     __zmlt__(matrix->me[i],scalar,out->me[i],(int)n);
  471.     /**************************************************
  472.       for ( j=0; j<n; j++ )
  473.       out->me[i][j] = scalar*matrix->me[i][j];
  474.       **************************************************/
  475.     return (out);
  476. }
  477.  
  478. /* zvm_mlt -- vector adjoint-matrix multiplication */
  479. ZVEC    *zvm_mlt(A,b,out)
  480. ZMAT    *A;
  481. ZVEC    *b,*out;
  482. {
  483.     u_int    j,m,n;
  484.     /* complex    sum,**A_v,*b_v; */
  485.     
  486.     if ( A==ZMNULL || b==ZVNULL )
  487.     error(E_NULL,"zvm_mlt");
  488.     if ( A->m != b->dim )
  489.     error(E_SIZES,"zvm_mlt");
  490.     if ( b == out )
  491.     error(E_INSITU,"zvm_mlt");
  492.     if ( out == ZVNULL || out->dim != A->n )
  493.     out = zv_resize(out,A->n);
  494.     
  495.     m = A->m;        n = A->n;
  496.     
  497.     zv_zero(out);
  498.     for ( j = 0; j < m; j++ )
  499.     if ( b->ve[j].re != 0.0 || b->ve[j].im != 0.0  )
  500.         __zmltadd__(out->ve,A->me[j],b->ve[j],(int)n,Z_CONJ);
  501.     /**************************************************
  502.       A_v = A->me;        b_v = b->ve;
  503.       for ( j=0; j<n; j++ )
  504.       {
  505.       sum = 0.0;
  506.       for ( i=0; i<m; i++ )
  507.       sum += b_v[