home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / memory < prev    next >
Encoding:
Text File  |  1994-04-04  |  19.5 KB  |  831 lines
  1.  & Zbigniew Leyk, all rights reserved.
  2. **
  3. **                 Meschach Library
  4. ** 
  5. ** This Meschach Library is provided "as is" without any express 
  6. ** or implied warranty of any kind with respect to this software. 
  7. ** In particular the authors shall not be liable for any direct, 
  8. ** indirect, special, incidental or consequential damages arising 
  9. ** in any way from use of the software.
  10. ** 
  11. ** Everyone is granted permission to copy, modify and redistribute this
  12. ** Meschach Library, provided:
  13. **  1.  All copies contain this copyright notice.
  14. **  2.  All modified copies shall carry a notice stating who
  15. **      made the last modification and the date of such modification.
  16. **  3.  No charge is made for this software or works derived from it.  
  17. **      This clause shall not be construed as constraining other software
  18. **      distributed on the same medium as this software, nor is a
  19. **      distribution fee considered a charge.
  20. **
  21. ***************************************************************************/
  22.  
  23.  
  24. /*
  25.     This file contains routines for import/exporting data to/from
  26.         MATLAB. The main routines are:
  27.             MAT *m_save(FILE *fp,MAT *A,char *name)
  28.             VEC *v_save(FILE *fp,VEC *x,char *name)
  29.             MAT *m_load(FILE *fp,char **name)
  30. */
  31.  
  32. #include        <stdio.h>
  33. #include        "matrix.h"
  34. #include    "matlab.h"
  35.  
  36. static char rcsid[] = "$Id: matlab.c,v 1.7 1994/01/13 05:30:09 des Exp $";
  37.  
  38. /* m_save -- save matrix in ".mat" file for MATLAB
  39.     -- returns matrix to be saved */
  40. MAT     *m_save(fp,A,name)
  41. FILE    *fp;
  42. MAT     *A;
  43. char    *name;
  44. {
  45.     int     i;
  46.     matlab  mat;
  47.  
  48.     if ( ! A )
  49.         error(E_NULL,"m_save");
  50.  
  51.     mat.type = 1000*MACH_ID + 100*ORDER + 10*PRECISION + 0;
  52.     mat.m = A->m;
  53.     mat.n = A->n;
  54.     mat.imag = FALSE;
  55.     mat.namlen = (name == (char *)NULL) ? 1 : strlen(name)+1;
  56.  
  57.     /* write header */
  58.     fwrite(&mat,sizeof(matlab),1,fp);
  59.     /* write name */
  60.     if ( name == (char *)NULL )
  61.         fwrite("",sizeof(char),1,fp);
  62.     else
  63.         fwrite(name,sizeof(char),(int)(mat.namlen),fp);
  64.     /* write actual data */
  65.     for ( i = 0; i < A->m; i++ )
  66.         fwrite(A->me[i],sizeof(Real),(int)(A->n),fp);
  67.  
  68.     return A;
  69. }
  70.  
  71.  
  72. /* v_save -- save vector in ".mat" file for MATLAB
  73.     -- saves it as a row vector
  74.     -- returns vector to be saved */
  75. VEC     *v_save(fp,x,name)
  76. FILE    *fp;
  77. VEC     *x;
  78. char    *name;
  79. {
  80.     matlab  mat;
  81.  
  82.     if ( ! x )
  83.         error(E_NULL,"v_save");
  84.  
  85.     mat.type = 1000*MACH_ID + 100*ORDER + 10*PRECISION + 0;
  86.     mat.m = x->dim;
  87.     mat.n = 1;
  88.     mat.imag = FALSE;
  89.     mat.namlen = (name == (char *)NULL) ? 1 : strlen(name)+1;
  90.  
  91.     /* write header */
  92.     fwrite(&mat,sizeof(matlab),1,fp);
  93.     /* write name */
  94.     if ( name == (char *)NULL )
  95.         fwrite("",sizeof(char),1,fp);
  96.     else
  97.         fwrite(name,sizeof(char),(int)(mat.namlen),fp);
  98.     /* write actual data */
  99.     fwrite(x->ve,sizeof(Real),(int)(x->dim),fp);
  100.  
  101.     return x;
  102. }
  103.  
  104. /* d_save -- save double in ".mat" file for MATLAB
  105.     -- saves it as a row vector
  106.     -- returns vector to be saved */
  107. double    d_save(fp,x,name)
  108. FILE    *fp;
  109. double    x;
  110. char    *name;
  111. {
  112.     matlab  mat;
  113.     Real x1 = x;
  114.  
  115.     mat.type = 1000*MACH_ID + 100*ORDER + 10*PRECISION + 0;
  116.     mat.m = 1;
  117.     mat.n = 1;
  118.     mat.imag = FALSE;
  119.     mat.namlen = (name == (char *)NULL) ? 1 : strlen(name)+1;
  120.  
  121.     /* write header */
  122.     fwrite(&mat,sizeof(matlab),1,fp);
  123.     /* write name */
  124.     if ( name == (char *)NULL )
  125.         fwrite("",sizeof(char),1,fp);
  126.     else
  127.         fwrite(name,sizeof(char),(int)(mat.namlen),fp);
  128.     /* write actual data */
  129.     fwrite(&x1,sizeof(Real),1,fp);
  130.  
  131.     return x;
  132. }
  133.  
  134. /* m_load -- loads in a ".mat" file variable as produced by MATLAB
  135.     -- matrix returned; imaginary parts ignored */
  136. MAT     *m_load(fp,name)
  137. FILE    *fp;
  138. char    **name;
  139. {
  140.     MAT     *A;
  141.     int     i;
  142.     int     m_flag, o_flag, p_flag, t_flag;
  143.     float   f_temp;
  144.     Real    d_temp;
  145.     matlab  mat;
  146.  
  147.     if ( fread(&mat,sizeof(matlab),1,fp) != 1 )
  148.         error(E_FORMAT,"m_load");
  149.     if ( mat.type >= 10000 )    /* don't load a sparse matrix! */
  150.         error(E_FORMAT,"m_load");
  151.     m_flag = (mat.type/1000) % 10;
  152.     o_flag = (mat.type/100) % 10;
  153.     p_flag = (mat.type/10) % 10;
  154.     t_flag = (mat.type) % 10;
  155.     if ( m_flag != MACH_ID )
  156.         error(E_FORMAT,"m_load");
  157.     if ( t_flag != 0 )
  158.         error(E_FORMAT,"m_load");
  159.     if ( p_flag != DOUBLE_PREC && p_flag != SINGLE_PREC )
  160.         error(E_FORMAT,"m_load");
  161.     *name = (char *)malloc((unsigned)(mat.namlen)+1);
  162.     if ( fread(*name,sizeof(char),(unsigned)(mat.namlen),fp) == 0 )
  163.         error(E_FORMAT,"m_load");
  164.     A = m_get((unsigned)(mat.m),(unsigned)(mat.n));
  165.     for ( i = 0; i < A->m*A->n; i++ )
  166.     {
  167.         if ( p_flag == DOUBLE_PREC )
  168.             fread(&d_temp,sizeof(double),1,fp);
  169.         else
  170.         {
  171.             fread(&f_temp,sizeof(float),1,fp);
  172.             d_temp = f_temp;
  173.         }
  174.         if ( o_flag == ROW_ORDER )
  175.             A->me[i / A->n][i % A->n] = d_temp;
  176.         else if ( o_flag == COL_ORDER )
  177.             A->me[i % A->m][i / A->m] = d_temp;
  178.         else
  179.             error(E_FORMAT,"m_load");
  180.     }
  181.  
  182.     if ( mat.imag )         /* skip imaginary part */
  183.     for ( i = 0; i < A->m*A->n; i++ )
  184.     {
  185.         if ( p_flag == DOUBLE_PREC )
  186.             fread(&d_temp,sizeof(double),1,fp);
  187.         else
  188.             fread(&f_temp,sizeof(float),1,fp);
  189.     }
  190.  
  191.     return A;
  192. }
  193.  
  194. FileDataŵmatopî.Eÿÿÿd…?◰.
  195. /**************************************************************************
  196. **
  197. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  198. **
  199. **                 Meschach Library
  200. ** 
  201. ** This Meschach Library is provided "as is" without any express 
  202. ** or implied warranty of any kind with respect to this software. 
  203. ** In particular the authors shall not be liable for any direct, 
  204. ** indirect, special, incidental or consequential damages arising 
  205. ** in any way from use of the software.
  206. ** 
  207. ** Everyone is granted permission to copy, modify and redistribute this
  208. ** Meschach Library, provided:
  209. **  1.  All copies contain this copyright notice.
  210. **  2.  All modified copies shall carry a notice stating who
  211. **      made the last modification and the date of such modification.
  212. **  3.  No charge is made for this software or works derived from it.  
  213. **      This clause shall not be construed as constraining other software
  214. **      distributed on the same medium as this software, nor is a
  215. **      distribution fee considered a charge.
  216. **
  217. ***************************************************************************/
  218.  
  219.  
  220. /* matop.c 1.3 11/25/87 */
  221.  
  222.  
  223. #include    <stdio.h>
  224. #include    "matrix.h"
  225.  
  226. static    char    rcsid[] = "$Id: matop.c,v 1.3 1994/01/13 05:30:28 des Exp $";
  227.  
  228.  
  229. /* m_add -- matrix addition -- may be in-situ */
  230. MAT    *m_add(mat1,mat2,out)
  231. MAT    *mat1,*mat2,*out;
  232. {
  233.     u_int    m,n,i;
  234.  
  235.     if ( mat1==(MAT *)NULL || mat2==(MAT *)NULL )
  236.         error(E_NULL,"m_add");
  237.     if ( mat1->m != mat2->m || mat1->n != mat2->n )
  238.         error(E_SIZES,"m_add");
  239.     if ( out==(MAT *)NULL || out->m != mat1->m || out->n != mat1->n )
  240.         out = m_resize(out,mat1->m,mat1->n);
  241.     m = mat1->m;    n = mat1->n;
  242.     for ( i=0; i<m; i++ )
  243.     {
  244.         __add__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  245.         /**************************************************
  246.         for ( j=0; j<n; j++ )
  247.             out->me[i][j] = mat1->me[i][j]+mat2->me[i][j];
  248.         **************************************************/
  249.     }
  250.  
  251.     return (out);
  252. }
  253.  
  254. /* m_sub -- matrix subtraction -- may be in-situ */
  255. MAT    *m_sub(mat1,mat2,out)
  256. MAT    *mat1,*mat2,*out;
  257. {
  258.     u_int    m,n,i;
  259.  
  260.     if ( mat1==(MAT *)NULL || mat2==(MAT *)NULL )
  261.         error(E_NULL,"m_sub");
  262.     if ( mat1->m != mat2->m || mat1->n != mat2->n )
  263.     e for any direct, 
  264. ** indirect, special, incidental or consequential damages arising 
  265. ** in any way from use of the software.
  266. ** 
  267. ** Everyone is granted permission to copy, modify and redistribute this
  268. ** Meschach Library, provided:
  269. **  1.  All copies contain this copyright notice.
  270. **  2.  All modified copies shall carry a notice stating who
  271. **      made the last modification and the date of such modification.
  272. **  3.  No charge is made for this software or works derived from it.  
  273. **      This clause shall not be construed as constraining other software
  274. **      distributed on the same medium as this software, nor is a
  275. **      distribution fee considered a charge.
  276. **
  277. ***************************************************************************/
  278.  
  279.  
  280. /*
  281.         Files for matrix computations
  282.  
  283.     Householder transformation file. Contains routines for calculating
  284.     householder transformations, applying them to vectors and matrices
  285.     by both row & column.
  286. */
  287.  
  288. /* hsehldr.c 1.3 10/8/87 */
  289. static    char    rcsid[] = "$Id: hsehldr.c,v 1.2 1994/01/13 05:36:29 des Exp $";
  290.  
  291. #include    <stdio.h>
  292. #include    <math.h>
  293. #include    "matrix.h"
  294. #include        "matrix2.h"
  295.  
  296.  
  297. /* hhvec -- calulates Householder vector to eliminate all entries after the
  298.     i0 entry of the vector vec. It is returned as out. May be in-situ */
  299. VEC    *hhvec(vec,i0,beta,out,newval)
  300. VEC    *vec,*out;
  301. u_int    i0;
  302. Real    *beta,*newval;
  303. {
  304.     Real    norm;
  305.  
  306.     out = _v_copy(vec,out,i0);
  307.     norm = sqrt(_in_prod(out,out,i0));
  308.     if ( norm <= 0.0 )
  309.     {
  310.         *beta = 0.0;
  311.         return (out);
  312.     }
  313.     *beta = 1.0/(norm * (norm+fabs(out->ve[i0])));
  314.     if ( out->ve[i0] > 0.0 )
  315.         *newval = -norm;
  316.     else
  317.         *newval = norm;
  318.     out->ve[i0] -= *newval;
  319.  
  320.     return (out);
  321. }
  322.  
  323. /* hhtrvec -- apply Householder transformation to vector -- may be in-situ */
  324. VEC    *hhtrvec(hh,beta,i0,in,out)
  325. VEC    *hh,*in,*out;    /* hh = Householder vector */
  326. u_int    i0;
  327. double    beta;
  328. {
  329.     Real    scale;
  330.     /* u_int    i; */
  331.  
  332.     if ( hh==(VEC *)NULL || in==(VEC *)NULL )
  333.         error(E_NULL,"hhtrvec");
  334.     if ( in->dim != hh->dim )
  335.         error(E_SIZES,"hhtrvec");
  336.     if ( i0 > in->dim )
  337.         error(E_BOUNDS,"hhtrvec");
  338.  
  339.     scale = beta*_in_prod(hh,in,i0);
  340.     out = v_copy(in,out);
  341.     __mltadd__(&(out->ve[i0]),&(hh->ve[i0]),-scale,(int)(in->dim-i0));
  342.     /************************************************************
  343.     for ( i=i0; i<in->dim; i++ )
  344.         out->ve[i] = in->ve[i] - scale*hh->ve[i];
  345.     ************************************************************/
  346.  
  347.     return (out);
  348. }
  349.  
  350. /* hhtrrows -- transform a matrix by a Householder vector by rows
  351.     starting at row i0 from column j0 -- in-situ */
  352. MAT    *hhtrrows(M,i0,j0,hh,beta)
  353. MAT    *M;
  354. u_int    i0, j0;
  355. VEC    *hh;
  356. double    beta;
  357. {
  358.     Real    ip, scale;
  359.     int    i /*, j */;
  360.  
  361.     if ( M==(MAT *)NULL || hh==(VEC *)NULL )
  362.         error(E_NULL,"hhtrrows");
  363.     if ( M->n != hh->dim )
  364.         error(E_RANGE,"hhtrrows");
  365.     if ( i0 > M->m || j0 > M->n )
  366.         error(E_BOUNDS,"hhtrrows");
  367.  
  368.     if ( beta == 0.0 )    return (M);
  369.  
  370.     /* for each row ... */
  371.     for ( i = i0; i < M->m; i++ )
  372.     {    /* compute inner product */
  373.         ip = __ip__(&(M->me[i][j0]),&(hh->ve[j0]),(int)(M->n-j0));
  374.         /**************************************************
  375.         ip = 0.0;
  376.         for ( j = j0; j < M->n; j++ )
  377.             ip += M->me[i][j]*hh->ve[j];
  378.         **************************************************/
  379.         scale = beta*ip;
  380.         if ( scale == 0.0 )
  381.             continue;
  382.  
  383.         /* do operation */
  384.         __mltadd__(&(M->me[i][j0]),&(hh->ve[j0]),-scale,
  385.                             (int)(M->n-j0));
  386.         /**************************************************
  387.         for ( j = j0; j < M->n; j++ )
  388.             M->me[i][j] -= scale*hh->ve[j];
  389.         **************************************************/
  390.     }
  391.  
  392.     return (M);
  393. }
  394.  
  395.  
  396. /* hhtrcols -- transform a matrix by a Householder vector by columns
  397.     starting at row i0 from column j0 -- in-situ */
  398. MAT    *hhtrcols(M,i0,j0,hh,beta)
  399. MAT    *M;
  400. u_int    i0, j0;
  401. VEC    *hh;
  402. double    beta;
  403. {
  404.     /* Real    ip, scale; */
  405.     int    i /*, k */;
  406.     static    VEC    *w = VNULL;
  407.  
  408.     if ( M==(MAT *)NULL || hh==(VEC *)NULL )
  409.         error(E_NULL,"hhtrcols");
  410.     if ( M->m != hh->dim )
  411.         error(E_SIZES,"hhtrcols");
  412.     if ( i0 > M->m || j0 > M->n )
  413.         error(E_BOUNDS,"hhtrcols");
  414.  
  415.     if ( beta == 0.0 )    return (M);
  416.  
  417.     w = v_resize(w,M->n);
  418.     MEM_STAT_REG(w,TYPE_VEC);
  419.     v_zero(w);
  420.  
  421.     for ( i = i0; i < M->m; i++ )
  422.         if ( hh->ve[i] != 0.0 )
  423.         __mltadd__(&(w->ve[j0]),&(M->me[i][j0]),hh->ve[i],
  424.                             (int)(M->n-j0));
  425.     for ( i = i0; i < M->m; i++ )
  426.         if ( hh->ve[i] != 0.0 )
  427.         __mltadd__(&(M->me[i][j0]),&(w->ve[j0]),-beta*hh->ve[i],
  428.                             (int)(M->n-j0));
  429.     return (M);
  430. }
  431.  
  432. FileDataŵinitEÿÿÿè%@Ðë
  433. /**************************************************************************
  434. **
  435. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  436. **
  437. **                 Meschach Library
  438. ** 
  439. ** This Meschach Library is provided "as is" without any express 
  440. ** or implied warranty of any kind with respect to this software. 
  441. ** In particular the authors shall not be liable for any direct, 
  442. ** indirect, special, incidental or consequential damages arising 
  443. ** in any way from use of the software.
  444. ** 
  445. ** Everyone is granted permission to copy, modify and redistribute this
  446. ** Meschach Library, provided:
  447. **  1.  All copies contain this copyright notice.
  448. **  2.  All modified copies shall carry a notice stating who
  449. **      made the last modification and the date of such modification.
  450. **  3.  No charge is made for this software or works derived from it.  
  451. **      This clause shall not be construed as constraining other software
  452. **      distributed on the same medium as this software, nor is a
  453. **      distribution fee considered a charge.
  454. **
  455. ***************************************************************************/
  456.  
  457.  
  458. /*
  459.     This is a file of routines for zero-ing, and initialising
  460.     vectors, matrices and permutations.
  461.     This is to be included in the matrix.a library
  462. */
  463.  
  464. static    char    rcsid[] = "$Id: init.c,v 1.6 1994/01/13 05:36:58 des Exp $";
  465.  
  466. #include    <stdio.h>
  467. #include    "matrix.h"
  468.  
  469. /* v_zero -- zero the vector x */
  470. VEC    *v_zero(x)
  471. VEC    *x;
  472. {
  473.     if ( x == VNULL )
  474.         error(E_NULL,"v_zero");
  475.  
  476.     __zero__(x->ve,x->dim);
  477.     /* for ( i = 0; i < x->dim; i++ )
  478.         x->ve[i] = 0.0; */
  479.  
  480.     return x;
  481. }
  482.  
  483.  
  484. /* iv_zero -- zero the vector ix */
  485. IVEC    *iv_zero(ix)
  486. IVEC    *ix;
  487. {
  488.    int i;
  489.    
  490.    if ( ix == IVNULL )
  491.      error(E_NULL,"iv_zero");
  492.    
  493.    for ( i = 0; i < ix->dim; i++ )
  494.      ix->ive[i] = 0; 
  495.    
  496.    return ix;
  497. }
  498.  
  499.  
  500. /* m_zero -- zero the matrix A */
  501. MAT    *m_zero(A)
  502. MAT    *A;
  503. {
  504.     int    i, A_m, A_n;
  505.     Real    **A_me;
  506.  
  507.     if ( A == MNULL )
  508.         error(E_NULL,"m_zero");
  509.  
  510.     A_m = A->m;    A_n = A->n;    A_me = A->me;
  511.     for ( i = 0; i < A_m; i++ )
  512.         __zero__(A_me[i],A_n);
  513.         /* for ( j = 0; j < A_n; j++ )
  514.             A_me[i][j] = 0.0; */
  515.  
  516.     return A;
  517. }
  518.  
  519. /* mat_id -- set A to being closest to identity matrix as possible
  520.     -- i.e. A[i][j] == 1 if i == j and 0 otherwise */
  521. MAT    *m_ident(A)
  522. MAT    *A;
  523. {
  524.     int    i, size;
  525.  
  526.     if ( A == MNULL )
  527.         error(E_NULL,"m_ident");
  528.  
  529.     m_zero(A);
  530.     size = min(A->m,A->n);
  531.     for ( i = 0; i < size; i++ )
  532.         A->me[i][i] = 1.0;
  533.  
  534.     return A;
  535. }
  536.     
  537. /* px_ident -- set px to identity permutation */
  538. PERM    *px_ident(px)
  539. PERM    *px;
  540. {
  541.     int    i, px_size;
  542.     u_int    *px_pe;
  543.  
  544.     if ( px == PNULL )
  545.         error(E_NULL,"px_ident");
  546.  
  547.     px_size = px->size;    px_pe = px->pe;
  548.     for ( i = 0; i < px_size; i++ )
  549.         px_pe[i] = i;
  550.  
  551.     return px;
  552. }
  553.  
  554. /* Pseudo random number generator data structures */
  555. /* Knuth's lagged Fibonacci-based generator: See "Seminumerical Algorithms:
  556.    The Art of Computer Programming" sections 3.2-3.3 */
  557.  
  558. #ifdef ANSI_C
  559. #ifndef LONG_MAX
  560. #include    <limits.h>
  561. #endif
  562. #endif
  563.  
  564. #ifdef LONG_MAX
  565. #define MODULUS    LONG_MAX
  566. #else
  567. #define MODULUS    1000000000L    /* assuming long's at least 32 bits long */
  568. #endif
  569. #define MZ    0L
  570.  
  571. static long mrand_list[56];
  572. static int  started = FALSE;
  573. static int  inext = 0, inextp = 31;
  574.  
  575.  
  576. /* mrand -- pseudo-random number generator */
  577. #ifdef ANSI_C
  578. double mrand(void)
  579. #else
  580. double mrand()
  581. #endif
  582. {
  583.     long    lval;
  584.     static Real  factor = 1.0/((Real)MODULUS);
  585.     
  586.     if ( ! started )
  587.     smrand(3127);
  588.     
  589.     inext = (inext >= 54) ? 0 : inext+1;
  590.     inextp = (inextp >= 54) ? 0 : inextp+1;
  591.  
  592.     lval = mrand_list[inext]-mrand_list[inextp];
  593.     if ( lval < 0L )
  594.     lval += MODULUS;
  595.     mrand_list[inext] = lval;
  596.     
  597.     return (double)lval*factor;
  598. }
  599.  
  600. /* mrandlist -- fills the array a[] with len random numbers */
  601. void    mrandlist(a, len)
  602. Real    a[];
  603. int    len;
  604. {
  605.     int        i;
  606.     long    lval;
  607.     static Real  factor = 1.0/((Real)MODULUS);
  608.     
  609.     if ( ! started )
  610.     smrand(3127);
  611.     
  612.     for ( i = 0; i < len; i++ )
  613.     {
  614.     inext = (inext >= 54) ? 0 : inext+1;
  615.     inextp = (inextp >= 54) ? 0 : inextp+1;
  616.     
  617.     lval = mrand_list[inext]-mrand_list[inextp];
  618.     if ( lval < 0L )
  619.         lval += MODULUS;
  620.     mrand_list[inext] = lval;
  621.     
  622.     a[i] = (Real)lval*factor;
  623.     }
  624. }
  625.  
  626.  
  627. /* smrand -- set seed for mrand() */
  628. void smrand(seed)
  629. int    seed;
  630. {
  631.     int        i;
  632.  
  633.     mrand_list[0] = (123413*seed) % MODULUS;
  634.     for ( i = 1; i < 55; i++ )
  635.     mrand_list[i] = (123413*mrand_list[i-1]) % MODULUS;
  636.  
  637.     started = TRUE;
  638.  
  639.     /* run mrand() through the list sufficient times to
  640.        thoroughly randomise the array */
  641.     for ( i = 0; i < 55*55; i++ )
  642.     mrand();
  643. }
  644. #undef MODULUS
  645. #undef MZ
  646. #undef FAC
  647.  
  648. /* v_rand -- initialises x to be a random vector, components
  649.     independently & uniformly ditributed between 0 and 1 */
  650. VEC    *v_rand(x)
  651. VEC    *x;
  652. {
  653.     /* int    i; */
  654.  
  655.     if ( ! x )
  656.         error(E_NULL,"v_rand");
  657.  
  658.     /* for ( i = 0; i < x->dim; i++ ) */
  659.         /* x->ve[i] = rand()/((Real)MAX_RAND); */
  660.         /* x->ve[i] = mrand(); */
  661.     mrandlist(x->ve,x->dim);
  662.  
  663.     return x;
  664. }
  665.  
  666. /* m_rand -- initialises A to be a random vector, components
  667.     independently & uniformly distributed between 0 and 1 */
  668. MAT    *m_rand(A)
  669. MAT    *A;
  670. {
  671.     int    i /* , j */;
  672.  
  673.     if ( ! A )
  674.         error(E_NULL,"m_rand");
  675.  
  676.     for ( i = 0; i < A->m; i++ )
  677.         /* for ( j = 0; j < A->n; j++ ) */
  678.             /* A->me[i][j] = rand()/((Real)MAX_RAND); */
  679.             /* A->me[i][j] = mrand(); */
  680.         mrandlist(A->me[i],A->n);
  681.  
  682.     return A;
  683. }
  684.  
  685. /* v_ones -- fills x with one's */
  686. VEC    *v_ones(x)
  687. VEC    *x;
  688. {
  689.     int    i;
  690.  
  691.     if ( ! x )
  692.         error(E_NULL,"v_ones");
  693.  
  694.     for ( i = 0; i < x->dim; i++ )
  695.         x->ve[i] = 1.0;
  696.  
  697.     return x;
  698. }
  699.  
  700. /* m_ones -- fills matrix with one's */
  701. MAT    *m_ones(A)
  702. MAT    *A;
  703. {
  704.     int    i, j;
  705.  
  706.     if ( ! A )
  707.         error(E_NULL,"m_ones");
  708.  
  709.     for ( i = 0; i < A->m; i++ )
  710.         for ( j = 0; j < A->n; j++ )
  711.             A->me[i][j] = 1.0;
  712.  
  713.     return A;
  714. }
  715.  
  716. /* v_count -- initialises x so that x->ve[i] == i */
  717. VEC    *v_count(x)
  718. VEC    *x;
  719. {
  720.     int    i;
  721.  
  722.     if ( ! x )
  723.         error(E_NULL,"v_count");
  724.  
  725.     for ( i = 0; i < x->d printf("# det_mant1=%g, det_expt1=%d\n",
  726.         det_mant1,det_expt1); */
  727.      /* printf("# det_mant2=%g, det_expt2=%d\n",
  728.         det_mant2,det_expt2); */
  729.      if ( det_mant1 == 0.0 )
  730.      {   /* multiple e-val of T */
  731.         err_est->ve[i] = 0.0;
  732.         continue;
  733.      }
  734.      else if ( det_mant2 == 0.0 )
  735.      {
  736.         err_est->ve[i] = HUGE;
  737.         continue;
  738.      }
  739.      if ( (det_expt1 + det_expt2) % 2 )
  740.        /* if odd... */
  741.        det_mant = sqrt(2.0*fabs(det_mant1*det_mant2));
  742.      else /* if even... */
  743.        det_mant = sqrt(fabs(det_mant1*det_mant2));
  744.      det_expt = (det_expt1+det_expt2)/2;
  745.      err_est->ve[i] = fabs(beta*
  746.                    ldexp(pb_mant/det_mant,pb_expt-det_expt));
  747.       }
  748.    }
  749.    
  750.    return a;
  751. }
  752.  
  753. /* iter_splanczos2 -- version of iter_lanczos2() that uses sparse matrix data
  754.    structure */
  755.  
  756. VEC    *iter_splanczos2(A,m,x0,evals,err_est)
  757. SPMAT    *A;
  758. int     m;
  759. VEC    *x0;        /* initial vector */
  760. VEC    *evals;        /* eigenvalue vector */
  761. VEC    *err_est;    /* error estimates of eigenvalues */
  762. {    
  763.    ITER *ip;
  764.    VEC *a;
  765.    
  766.    ip = iter_get(0,0);
  767.    ip->Ax = (Fun_Ax) sp_mv_mlt;
  768.    ip->A_par = (void *) A;
  769.    ip->x = x0;
  770.    ip->k = m;
  771.    a = iter_lanczos2(ip,evals,err_est);    
  772.    ip->shared_x = ip->shared_b = TRUE;
  773.    iter_free(ip);   /* release only ITER structure */
  774.    return a;
  775. }
  776.  
  777.  
  778.  
  779.  
  780. /*
  781.   Conjugate gradient method
  782.   Another variant - mainly for testing
  783.   */
  784.  
  785. VEC  *iter_cg1(ip)
  786. ITER *ip;
  787. {
  788.    static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
  789.    Real    alpha;
  790.    double inner,nres;
  791.    VEC *rr;   /* rr == r or rr == z */
  792.    
  793.    if (ip == INULL)
  794.      error(E_NULL,"iter_cg");
  795.    if (!ip->Ax || !ip->b)
  796.      error(E_NULL,"iter_cg");
  797.    if ( ip->x == ip->b )
  798.      error(E_INSITU,"iter_cg");
  799.    if (!ip->stop_crit)
  800.      error(E_NULL,"iter_cg");
  801.    
  802.    if ( ip->eps <= 0.0 )
  803.      ip->eps = MACHEPS;
  804.    
  805.    r = v_resize(r,ip->b->dim);
  806.    p = v_resize(p,ip->b->dim);
  807.    q = v_resize(q,ip->b->dim);
  808.    
  809.    MEM_STAT_REG(r,TYPE_VEC);
  810.    MEM_STAT_REG(p,TYPE_VEC);
  811.    MEM_STAT_REG(q,TYPE_VEC);
  812.    
  813.    if (ip->Bx != (Fun_Ax)NULL) {
  814.       z = v_resize(z,ip->b->dim);
  815.       MEM_STAT_REG(z,TYPE_VEC    error(E_SIZES,"m_sub");
  816.     if ( out==(MAT *)NULL || out->m != mat1->m || out->n != mat1->n )
  817.         out = m_resize(out,mat1->m,mat1->n);
  818.     m = mat1->m;    n = mat1->n;
  819.     for ( i=0; i<m; i++ )
  820.     {
  821.         __sub__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  822.         /**************************************************
  823.         for ( j=0; j<n; j++ )
  824.             out->me[i][j] = mat1->me[i][j]-mat2->me[i][j];
  825.         **************************************************/
  826.     }
  827.  
  828.     return (out);
  829. }
  830.  
  831. /* m_mlt -- matrix-m