Back to index

plt-scheme  4.2.1
Classes | Defines | Typedefs | Functions | Variables
matrix.c File Reference
#include "fit.h"
#include "matrix.h"
#include <math.h>

Go to the source code of this file.

Classes

struct  Node

Defines

#define NULL   0
#define null   0
#define Swap(a, b)   {double temp = (a); (a) = (b); (b) = temp;}
#define WINZIG   1e-30

Typedefs

typedef struct Node Node

Functions

voidmy_gc_malloc (int size)
void gc_cleanup ()
static int fsign (double x)
MZ_DLLEXPORT double * vec (int n)
MZ_DLLEXPORT double ** matr (int rows, int cols)
void free_matr (double **m)
MZ_DLLEXPORT double * redim_vec (double **v, int n)
MZ_DLLEXPORT void redim_ivec (int **v, int n)
MZ_DLLEXPORT void Givens (double **C, double *d, double *x, double *r, int N, int n, int want_r)
MZ_DLLEXPORT void Invert_RtR (double **R, double **I, int n)

Variables

static char * RCSid = "$Id: matrix.c,v 1.4 2005/03/15 23:21:26 eli Exp $"
Nodehead = null

Class Documentation

struct Node

Definition at line 46 of file matrix.c.

Collaboration diagram for Node:
Class Members
struct Node * next
void * ptr

Define Documentation

#define NULL   0

Definition at line 38 of file matrix.c.

#define null   0

Definition at line 39 of file matrix.c.

#define Swap (   a,
  b 
)    {double temp = (a); (a) = (b); (b) = temp;}

Definition at line 75 of file matrix.c.

#define WINZIG   1e-30

Definition at line 76 of file matrix.c.


Typedef Documentation

typedef struct Node Node

Function Documentation

void free_matr ( double **  m)

Definition at line 121 of file matrix.c.

{
    free (m[0]);
    free (m);
}
static int fsign ( double  x) [static]

Definition at line 159 of file matrix.c.

{
    return( x>0 ? 1 : (x < 0) ? -1 : 0) ;
}

Here is the caller graph for this function:

Definition at line 62 of file matrix.c.

                 {
  while(head) {
    Node * current = head;
    head = current->next;
    free(current->ptr);
    free(current);
  }
}

Here is the caller graph for this function:

MZ_DLLEXPORT void Givens ( double **  C,
double *  d,
double *  x,
double *  r,
int  N,
int  n,
int  want_r 
)

Definition at line 181 of file matrix.c.

{
    int i, j, k;
    double w, gamma, sigma, rho, temp;
    double epsilon = 1e-5; /* FIXME (?)*/

/* 
 * First, construct QR decomposition of C, by 'rotating away'
 * all elements of C below the diagonal. The rotations are
 * stored in place as Givens coefficients rho.
 * Vector d is also rotated in this same turn, if it exists 
 */
    for (j = 0; j<n; j++) 
       for (i = j+1; i<N; i++) 
           if (C[i][j]) {
              if (fabs(C[j][j])<epsilon*fabs(C[i][j])) { /* find the rotation parameters */
                  w = -C[i][j];
                  gamma = 0;
                  sigma = 1;
                  rho = 1;
              } else {
                  w = fsign(C[j][j])*sqrt(C[j][j]*C[j][j] + C[i][j]*C[i][j]);
                  if (w == 0) {
                    //                    Eex3 ( "w = 0 in Givens();  Cjj = %g,  Cij = %g", C[j][j], C[i][j]);
                  }
                  gamma = C[j][j]/w;
                  sigma = -C[i][j]/w;
                  rho = (fabs(sigma)<gamma) ? sigma : fsign(sigma)/gamma;
              }
              C[j][j] = w;
              C[i][j] = rho;           /* store rho in place, for later use */
              for (k = j+1; k<n; k++) {   /* rotation on index pair (i,j) */
                  temp =    gamma*C[j][k] - sigma*C[i][k];
                  C[i][k] = sigma*C[j][k] + gamma*C[i][k];
                  C[j][k] = temp;
                  
              }
              if (d) {               /* if no d vector given, don't use it */
                  temp = gamma*d[j] - sigma*d[i];  /* rotate d */
                  d[i] = sigma*d[j] + gamma*d[i];
                  d[j] = temp;
               }
           }
    if (!d)               /* stop here if no d was specified */
         return;

    for (i = n-1; i >= 0; i--) {   /* solve R*x+d = 0, by backsubstitution */
        double s = d[i];
        r[i] = 0;              /* ... and also set r[i] = 0 for i<n */
        for (k = i+1; k<n; k++) 
            s += C[i][k]*x[k];
       if (C[i][i] == 0) {
         //Eex ( "Singular matrix in Givens()");
       }
        x[i] = - s / C[i][i];
       }
    for (i = n; i < N; i++) 
       r[i] = d[i];         /* set the other r[i] to d[i] */
       
    if (!want_r)            /* if r isn't needed, stop here */
       return;
       
    /* rotate back the r vector */
    for (j = n-1; j >= 0; j--)
       for (i = N-1; i >= 0; i--) {
           if ((rho = C[i][j]) == 1) { /* reconstruct gamma, sigma from stored rho */
              gamma = 0;
              sigma = 1;
           } else if (fabs(rho)<1) {
              sigma = rho; 
              gamma = sqrt(1-sigma*sigma);
           } else {
              gamma = 1/fabs(rho);
              sigma = fsign(rho)*sqrt(1-gamma*gamma);
           }
           temp = gamma*r[j] + sigma*r[i];       /* rotate back indices (i,j) */
           r[i] = -sigma*r[j] + gamma*r[i];
           r[j] = temp;
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

MZ_DLLEXPORT void Invert_RtR ( double **  R,
double **  I,
int  n 
)

Definition at line 280 of file matrix.c.

{
  int i, j, k;

  /* fill in the I matrix, and check R for regularity : */

  for (i = 0; i<n; i++) {
    for (j = 0; j<i; j++)  /* upper triangle isn't needed */
      I[i][j] = 0;
    I[i][i] = 1;
    if (! R[i][i])
      {
      //      Eex ("Singular matrix in Invert_RtR");
      }
      }
  
  /* Forward substitution: Solve R^T * B = I, store B in place of I */
  
  for (k = 0; k<n; k++) 
    for (i = k; i<n; i++) {  /* upper half needn't be computed */
      double s = I[i][k];
      for (j = k; j<i; j++)  /* for j<k, I[j][k] always stays zero! */
       s -= R[j][i] * I[j][k];
      I[i][k] = s / R[i][i];
    }

  /* Backward substitution: Solve R * A = B, store A in place of B */

  for (k = 0; k<n; k++)
    for (i = n-1; i >= k; i--) {  /* don't compute upper triangle of A */
      double s = I[i][k];
      for (j = i+1; j<n; j++)
       s -= R[i][j] * I[j][k];
      I[i][k] = s / R[i][i]; 
    }
}

Here is the caller graph for this function:

MZ_DLLEXPORT double** matr ( int  rows,
int  cols 
)

Definition at line 102 of file matrix.c.

{
    /* allocates a double matrix */

    register int i;
    register double **m;

    if ( rows < 1  ||  cols < 1 )
        return NULL;
    m = (double **) my_gc_malloc (rows * sizeof(double *));
    m[0] = (double *) my_gc_malloc (rows * cols * sizeof(double));
    for ( i = 1; i<rows ; i++ )
       m[i] = m[i-1] + cols;
    return m;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void* my_gc_malloc ( int  size)

Definition at line 53 of file matrix.c.

                              {
  void * ptr = malloc(size);
  Node * n = (Node *)malloc(sizeof (Node));
  n->ptr = ptr;
  n->next = head;
  head = n;
  return ptr;
}

Here is the caller graph for this function:

MZ_DLLEXPORT void redim_ivec ( int **  v,
int  n 
)

Definition at line 142 of file matrix.c.

{
    if ( n < 1 ) {
       *v = NULL;
       return;
    }
    *v = (int *) my_gc_malloc ( n * sizeof(int));
}

Here is the call graph for this function:

MZ_DLLEXPORT double* redim_vec ( double **  v,
int  n 
)

Definition at line 130 of file matrix.c.

{
    if ( n < 1 ) 
      *v = NULL;
    else       
      *v = (double *) my_gc_malloc( n * sizeof(double));
    return *v;
}

Here is the call graph for this function:

MZ_DLLEXPORT double* vec ( int  n)

Definition at line 89 of file matrix.c.

{
    /* allocates a double vector with n elements */
    double *dp;
    if( n < 1 )
       return (double *) NULL;
    dp = (double *) my_gc_malloc (n * sizeof(double));
    return dp;
}

Here is the call graph for this function:


Variable Documentation

Definition at line 51 of file matrix.c.

char* RCSid = "$Id: matrix.c,v 1.4 2005/03/15 23:21:26 eli Exp $" [static]

Definition at line 2 of file matrix.c.