home *** CD-ROM | disk | FTP | other *** search
/ vis-ftp.cs.umass.edu / vis-ftp.cs.umass.edu.tar / vis-ftp.cs.umass.edu / pub / Software / ASCENDER / ascender.tar.Z / ascender.tar / Triangulate / array_utils.c next >
C/C++ Source or Header  |  1995-04-13  |  13KB  |  556 lines
  1. /* @(#)array_utils.c    1.15 12/20/94 */
  2.  
  3. /*>>
  4. Utility routines for the Numerical recipes cookbook programs 
  5. <<*/
  6.  
  7. #include "cvar.h"
  8. #include "rh_util.h"
  9. #include "alloc.h"
  10. #include <math.h>
  11. #include "newgaussj.s"
  12. #include "newgaussj2.s"
  13.  
  14. FUNCTION_DEF ( char **alloca_matrix, (
  15.         int nrl,
  16.         int nrh,
  17.         int ncl,
  18.         int nch,
  19.         int entry_size,
  20.         char **m))
  21.    {
  22.    /* Allocates a matrix with the range [nrl .. nrh][ncl .. nch] */
  23.    int i;
  24.    char *start;            /* Address of the [0][0] entry */
  25.  
  26.    int rdim = nrh - nrl + 1;    /* # entries in a row */
  27.    int cdim = nch - ncl + 1;    /* # entries in a column */
  28.  
  29.    /* Adjust the start of arrays, and adjust m */
  30.    start = (char *) (m + rdim);    /* Count forward rdim (MTYPE *) */
  31.    start -= (ncl + nrl * cdim) * entry_size; /* Count back units of entry_size*/
  32.    m -= nrl;
  33.  
  34.    /* Fill out the addresses */
  35.    for (i=nrl; i<=nrh; i++)
  36.       m[i] = &(start[i*cdim*entry_size]);
  37.  
  38.    /* Return pointer to the pointer array */
  39.    return (m);
  40.    }
  41.  
  42. FUNCTION_DEF ( char **malloc_matrix, (
  43.         int nrl,
  44.         int nrh,
  45.         int ncl,
  46.         int nch,
  47.         int entry_size))
  48.    {
  49.    /* Allocates a matrix with the range [nrl .. nrh][ncl .. nch] */
  50.    int i;
  51.    char **m;
  52.    char *start;            /* Address of the [0][0] entry */
  53.  
  54.    int rdim = nrh - nrl + 1;    /* # entries in a row */
  55.    int cdim = nch - ncl + 1;    /* # entries in a column */
  56.    int nptrs = rdim;
  57.    int size;
  58.  
  59.    /*
  60.     * We must make sure that double words are allocated on double word
  61.     * boundaries
  62.     */
  63.    /* We must make sure that we have an even number of pointers */
  64.    if (entry_size > sizeof(char *) && nptrs % 2 != 0)
  65.       nptrs += 1;
  66.       
  67.    /* Now, compute the size of space to be allocated */
  68.    size = nptrs*sizeof(char *) + rdim*cdim*entry_size;
  69.  
  70.    /* Allocate m */
  71.    m = (char **) MALLOC(size);
  72.    start = (char *) (m + nptrs);             /* Count forward nptrs (MTYPE *) */
  73.    start -= (ncl + nrl * cdim) * entry_size; /* Count back units of entry_size*/
  74.    m -= nrl;
  75.  
  76.    /* Fill out the addresses */
  77.    for (i=nrl; i<=nrh; i++)
  78.       m[i] = &(start[i*cdim*entry_size]);
  79.  
  80.    /* Return pointer to the pointer array */
  81.    return (m);
  82.    }
  83.  
  84. FUNCTION_DEF ( void fprint_intmatrix, (
  85.         char *format,
  86.         FILE *afile,
  87.         int **a,
  88.         int nrl,
  89.         int nrh,
  90.         int ncl,
  91.         int nch))
  92.    {
  93.    /* Print out a matrix of int values */
  94.    int i, j;
  95.    for (i=nrl; i<=nrh; i++)
  96.       {
  97.       for (j=ncl; j<=nch; j++)
  98.          fprintf (afile, format, a[i][j]);
  99.       fprintf (afile, "\n");
  100.       }
  101.    fprintf (afile, "\n");
  102.    }
  103.  
  104. FUNCTION_DEF ( void fprint_shortmatrix, (
  105.         char *format,
  106.         FILE *afile,
  107.         short **a,
  108.         int nrl,
  109.         int nrh,
  110.         int ncl,
  111.         int nch))
  112.    {
  113.    /* Print out a matrix of short values */
  114.    int i, j;
  115.    for (i=nrl; i<=nrh; i++)
  116.       {
  117.       for (j=ncl; j<=nch; j++)
  118.          fprintf (afile, format, a[i][j]);
  119.       fprintf (afile, "\n");
  120.       }
  121.    fprintf (afile, "\n");
  122.    }
  123.  
  124. FUNCTION_DEF ( void fprint_floatmatrix, (
  125.         char *format,
  126.         FILE *afile,
  127.         float **a,
  128.         int nrl,
  129.         int nrh,
  130.         int ncl,
  131.         int nch))
  132.    {
  133.    /* Print out a matrix of float values */
  134.    int i, j;
  135.    for (i=nrl; i<=nrh; i++)
  136.       {
  137.       for (j=ncl; j<=nch; j++)
  138.          fprintf (afile, format, a[i][j]);
  139.       fprintf (afile, "\n");
  140.       }
  141.    fprintf (afile, "\n");
  142.    }
  143.  
  144. FUNCTION_DEF ( void fprint_doublematrix, (
  145.         char *format,
  146.         FILE *afile,
  147.         double **a,
  148.         int nrl,
  149.         int nrh,
  150.         int ncl,
  151.         int nch))
  152.    {
  153.    /* Print out a matrix of double values */
  154.    int i, j;
  155.    for (i=nrl; i<=nrh; i++)
  156.       {
  157.       for (j=ncl; j<=nch; j++)
  158.          fprintf (afile, format, a[i][j]);
  159.       fprintf (afile, "\n");
  160.       }
  161.    fprintf (afile, "\n");
  162.    }
  163.  
  164. FUNCTION_DEF ( void fprint_matrix, (
  165.         char *format,
  166.         FILE *afile,
  167.         double **a,
  168.         int nrl,
  169.         int nrh,
  170.         int ncl,
  171.         int nch))
  172.    {
  173.    /* Print out a matrix of double values */
  174.    int i, j;
  175.    for (i=nrl; i<=nrh; i++)
  176.       {
  177.       for (j=ncl; j<=nch; j++)
  178.          fprintf (afile, format, a[i][j]);
  179.       fprintf (afile, "\n");
  180.       }
  181.    fprintf (afile, "\n");
  182.    }
  183.  
  184. FUNCTION_DEF ( void print_matrix, (
  185.         FILE *afile,
  186.         double **a,
  187.         int nrl,
  188.         int nrh,
  189.         int ncl,
  190.         int nch))
  191.    {
  192.    /* Print out a matrix of double values */
  193.    int i, j;
  194.    for (i=nrl; i<=nrh; i++)
  195.       {
  196.       int count = 0;
  197.       for (j=ncl; j<=nch; j++)
  198.      {
  199.          if (++count == 8)
  200.         {
  201.         fprintf (afile, "\n\t");
  202.         count = 1;
  203.         }
  204.          fprintf (afile,"%12.6f ", a[i][j]);
  205.      }
  206.       fprintf (afile, "\n");
  207.       }
  208.    fprintf (afile, "\n");
  209.    }
  210.  
  211. FUNCTION_DEF ( void fprint_vector, (
  212.         char *format,
  213.         FILE *afile,
  214.         double *v,
  215.         int nrl,
  216.         int nrh))
  217.    {
  218.    /* Print out a vector of double values */
  219.    int i;
  220.    for (i=nrl; i<=nrh; i++)
  221.       fprintf (afile, format, v[i]);
  222.    fprintf (afile, "\n");
  223.    }
  224.  
  225. FUNCTION_DEF ( void print_vector, (
  226.         FILE *afile,
  227.         double *v,
  228.         int nrl,
  229.         int nrh))
  230.    {
  231.    /* Print out a vector of double values */
  232.    int i;
  233.    int count = 0;
  234.    for (i=nrl; i<=nrh; i++)
  235.       {
  236.       if (++count == 8)
  237.      {
  238.      fprintf (afile, "\n\t");
  239.      count = 1;
  240.      }
  241.       fprintf (afile, "%9.6f ", v[i]);
  242.       }
  243.    fprintf (afile, "\n\n");
  244.    }
  245.  
  246. FUNCTION_DEF ( void matrix_prod, (
  247.         double **A,
  248.         double **B,
  249.         int l,
  250.         int m,
  251.         int n,
  252.         double **X))
  253.    {
  254.    /* Multiplies a k x m matrix (A) by an m x n matrix (B) to get matrix X */
  255.    /* All indices start at 1 */
  256.    /* The goal matrix may be the same as one of the operands */
  257.    int i, j, k;
  258.    double **T = MATRIX (1, l, 1, n, double);     /* Temporary product */
  259.  
  260.    /* Clear the matrix */
  261.    for (i=1; i<=l; i++)
  262.       for (j=1; j<=n; j++)
  263.          T[i][j] = 0.0;
  264.  
  265.    /* Do the multiplication */
  266.    for (i=1; i<=l; i++)
  267.       for (j=1; j<=n; j++)
  268.          for (k=1; k<=m; k++)
  269.             T[i][j] += A[i][k] * B[k][j];
  270.  
  271.    /* Clear the matrix */
  272.    for (i=1; i<=l; i++)
  273.       for (j=1; j<=n; j++)
  274.          X[i][j] = T[i][j];
  275.  
  276.    /* Free T */
  277.    FREE (T, 1);
  278.    }
  279.  
  280.  
  281. FUNCTION_DEF ( void vector_matrix_prod, (
  282.         double *b,
  283.         double **A,
  284.         int l,
  285.         int m,
  286.         double *x))
  287.    {
  288.    /* Multiplies a vector (b) by an  l x m matrix (A) to get vector x */
  289.    /* All indices start at 1 */
  290.    /* The goal vector may be the same as b */
  291.    int i, j, k;
  292.    double *T = VECTOR (1, m, double);     /* Temporary product */
  293.  
  294.    /* Clear the vector */
  295.    for (i=1; i<=m; i++)
  296.       T[i] = 0.0;
  297.  
  298.    /* Do the multiplication */
  299.    for (i=1; i<=m; i++)
  300.       for (k=1; k<=l; k++)
  301.          T[i] += b[k] * A[k][i];
  302.  
  303.    /* Transfer the vector */
  304.    for (i=1; i<=m; i++)
  305.       x[i] = T[i];
  306.  
  307.    /* Free T */
  308.    FREE (T, 1);
  309.    }
  310.  
  311. FUNCTION_DEF ( void matrix_vector_prod, (
  312.         double **A,
  313.         double *b,
  314.         int l,
  315.         int m,
  316.         double *x))
  317.    {
  318.    /* Multiplies a l x m matrix (A) by an m vector (b) to get vector x */
  319.    /* All indices start at 1 */
  320.    /* The goal vector may be the same as b */
  321.    int i, j, k;
  322.    double *T = VECTOR (1, l, double);     /* Temporary product */
  323.  
  324.    /* Clear the vector */
  325.    for (i=1; i<=l; i++)
  326.       T[i] = 0.0;
  327.  
  328.    /* Do the multiplication */
  329.    for (i=1; i<=l; i++)
  330.       for (k=1; k<=m; k++)
  331.          T[i] += A[i][k] * b[k];
  332.  
  333.    /* Transfer the vector */
  334.    for (i=1; i<=l; i++)
  335.       x[i] = T[i];
  336.  
  337.    /* Free T */
  338.    FREE (T, 1);
  339.    }
  340.  
  341.  
  342. FUNCTION_DEF ( double inner_product, (
  343.         double *A,
  344.         double *B,
  345.         int n))
  346.    {
  347.    /* Returns the inner product of two vectors */
  348.    int i;
  349.    double prod = 0.0;
  350.  
  351.    /* Do the multiplication */
  352.    for (i=1; i<=n; i++)
  353.       prod += A[i] * B[i];
  354.  
  355.    /* Return the product */
  356.    return (prod);
  357.    }
  358.  
  359. FUNCTION_DEF ( double vector_length, (
  360.         double *A,
  361.         int n))
  362.    {
  363.    /* Returns the length of the vector */
  364.    return (sqrt (inner_product (A, A, n)));
  365.    }
  366.  
  367. FUNCTION_DEF ( void transpose, (
  368.    double **A,        /* The matrix to be transposed */
  369.    double **At,        /* Its transpose (returned) */
  370.    int nrl, int nrh,    /* X dimension */
  371.    int ncl, int nch    /* Y dimension */
  372.    ))
  373.    {
  374.    /* Does the transpose of a matrix.          */
  375.    /* Goal and origin may be the same.        */
  376.    int i, j;
  377.    double **temp;
  378.    int del = 0;
  379.    
  380.    /* First does the diagonal entries */
  381.    if (A == At)
  382.       {
  383.       temp = MATRIX (nrl, nrh, ncl, nch, double);
  384.       del = 1;
  385.       for (i=nrl; i<=nrh; i++) for (j=1; j<=nch; j++)
  386.       temp[i][j] = A[i][j];
  387.       }
  388.    else temp = A;
  389.  
  390.    /* Now swaps the off-diagonal entries */
  391.    for (i=nrl; i<=nrh; i++) for (j=ncl; j<=nch; j++)
  392.       At[j][i] = temp[i][j];
  393.  
  394.    /* Free the matrix if required */
  395.    if(del) FREE(temp, nrl);
  396.    }
  397.  
  398.  
  399. FUNCTION_DEF ( int matrix_inverse, (
  400.    double **A,        /* The matrix to be inverted */
  401.    double **Ainv,    /* The inverse computed */
  402.    int n        /* Dimension of the matrix */
  403.    ))
  404.    {
  405.    /* Inverts the matrix */
  406.    /* Returns 0 if matrix is not invertible */
  407.    int i, j;
  408.    
  409.    /* First copy the matrix to Ainv */
  410.    for (i=1; i<=n; i++) for (j=1; j<=n; j++)
  411.       Ainv[i][j] = A[i][j];
  412.  
  413.    /* Now invert it on the spot */
  414.    return gaussj2_check (Ainv, n, (double **)0, 0);
  415.    }
  416.  
  417. FUNCTION_DEF ( char **sub_matrix, (
  418.    char **A,            /* Base array */
  419.    int nrl, int ncl,        /* Origin of the new array */
  420.    int rdim,            /* Number of rows of the sub-array */
  421.    int orgi, int orgj,        /* Starting index for new array */
  422.    int entry_size        /* Size of an entry */
  423.    ))
  424.    {
  425.    /* Allocates a sub-matrix of A with the given range */
  426.    int i;
  427.    char **m;
  428.    int offset;
  429.  
  430.    int size = rdim*sizeof(char *);
  431.  
  432.    /* Allocate m */
  433.    m = (char **) MALLOC(size);
  434.    m -= nrl ;
  435.  
  436.    /* Fill out the addresses */
  437.    offset = (ncl-orgj) * (entry_size/sizeof(char));
  438.    for (i=nrl; i<nrl+rdim; i++)
  439.       m[i] = A[i] + offset;
  440.  
  441.    /* Return pointer to the pointer array */
  442.    return (m+nrl-orgi);
  443.    }
  444.  
  445. FUNCTION_DEF ( double det3x3, (
  446.         double **A))
  447.    {
  448.    /* Takes a 3 by 3 determinant */
  449.    return A[1][1]*A[2][2]*A[3][3] + A[1][2]*A[2][3]*A[3][1] + 
  450.       A[1][3]*A[2][1]*A[3][2] - A[1][1]*A[2][3]*A[3][2] - 
  451.       A[1][2]*A[2][1]*A[3][3] - A[1][3]*A[2][2]*A[3][1];
  452.    }
  453.  
  454. FUNCTION_DEF ( void cofactor_matrix_3x3, (
  455.         double **A,
  456.         double **Aadj))
  457.    {
  458.    /* Takes the cofactor matrix of a given 3x3 matrix */
  459.    int i, j;
  460.    
  461.    /* First, take a copy of A */
  462.    double **A2 = MATRIX (1, 3, 1, 3, double);
  463.    for (i=1; i<=3; i++) for (j=1; j<=3; j++)
  464.       A2[i][j] = A[i][j];
  465.  
  466.    /* Now take the cofactors */
  467.    Aadj[1][1] = A2[2][2]*A2[3][3] - A2[2][3]*A2[3][2];
  468.    Aadj[1][2] = A2[2][3]*A2[3][1] - A2[2][1]*A2[3][3];
  469.    Aadj[1][3] = A2[2][1]*A2[3][2] - A2[2][2]*A2[3][1];
  470.  
  471.    Aadj[2][1] = A2[3][2]*A2[1][3] - A2[3][3]*A2[1][2];
  472.    Aadj[2][2] = A2[3][3]*A2[1][1] - A2[3][1]*A2[1][3];
  473.    Aadj[2][3] = A2[3][1]*A2[1][2] - A2[3][2]*A2[1][1];
  474.  
  475.    Aadj[3][1] = A2[1][2]*A2[2][3] - A2[1][3]*A2[2][2];
  476.    Aadj[3][2] = A2[1][3]*A2[2][1] - A2[1][1]*A2[2][3];
  477.    Aadj[3][3] = A2[1][1]*A2[2][2] - A2[1][2]*A2[2][1];
  478.  
  479.    /* Free the temporary matrix */
  480.    FREE (A2, 1);
  481.    }
  482.  
  483. FUNCTION_DEF ( void matrix_copy, (
  484.    double **A, double **B,    /* The matrices to copy from and to */
  485.    int i0, int i1,        /* Row bounds */
  486.    int j0, int j1        /* Column bounds */
  487.    ))
  488.    {
  489.    /* Copies a matrix */
  490.    int i, j;
  491.    for (i=i0; i<=i1; i++)  
  492.       for (j=j0; j<=j1; j++)
  493.      B[i][j] = A[i][j];
  494.    }
  495.    
  496. FUNCTION_DEF ( void cross_product, (
  497.    double *A, double *B,    /* Two vectors 1..3 */
  498.    double *C            /* Their cross product */
  499.    ))
  500.    {
  501.    C[1] = A[2]*B[3] - A[3]*B[2];
  502.    C[2] = A[3]*B[1] - A[1]*B[3];
  503.    C[3] = A[1]*B[2] - A[2]*B[1];
  504.    }
  505.  
  506. FUNCTION_DEF ( int lin_solve, (
  507.    double **A,        /* The matrix of coefficients */
  508.    double *x,        /* The vector of unknowns */
  509.    double *b,        /* The goal vector */
  510.    int n        /* Size of the system */
  511.    ))
  512.    {
  513.    /* Solves a square linear system of equations */
  514.    int i, j;
  515.    int return_val;
  516.    double **A2 = MATRIX (1, n, 1, n, double);
  517.    double *b2  = VECTOR (1, n, double);
  518.  
  519.    /* Copy the values */
  520.    for (i=1; i<=n; i++) for (j=1; j<=n; j++)
  521.       A2[i][j] = A[i][j];
  522.    for (i=1; i<=n; i++)
  523.       b2[i] = b[i];
  524.  
  525.    /* Now solve using gaussj */
  526.    return_val = gaussj_check (A2, n, b2);
  527.  
  528.    /* Now copy the values */
  529.    if (return_val)
  530.       for (i=1; i<=n; i++)
  531.          x[i] = b2[i];
  532.  
  533.    /* Free the temporary matrices */
  534.    FREE (A2, 1);
  535.    FREE (b2, 1);
  536.  
  537.    /* Return the appropriate value */
  538.    return return_val;
  539.    }
  540.    
  541. FUNCTION_DEF ( void identity_matrix, (
  542.         double **A,
  543.         int dim))
  544.    { 
  545.    /* Initializes the matrix A to be the identity */
  546.    int i, j;
  547.    for (i=1; i<=dim; i++) 
  548.       { 
  549.       A[i][i] = 1.0; 
  550.       for (j=i+1; j<=dim; j++)
  551.          A[i][j] = A[j][i] = 0.0;
  552.       } 
  553.    } 
  554.  
  555.  
  556.