Code owners
Assign users and groups as approvers for specific file changes. Learn more.
nrcheck.c 22.59 KiB
/* nrcheck.c **************************************************************
NR Library Package - Wrapper functions for thorough checking
James Trevelyan, University of Western Australia
Revision 2 January 1996
To use this to add safty checks to your program, add the
module NRCHECK.C to your project, and #define NR_CHECK before
including "nrlin.h" in your source file.
**************************************************************************/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "nrlin.h"
/* #define TEST */ /* use to activate test program at end of source */
#define valid_dvector(v,nl,nh) ((v!=(double *)NULL) && \
((v[(nl)-1]==-322.0)&&(v[(nh)+1]==-722.0)))
#define valid_vector(v,nl,nh) ((v!=(float *)NULL) && \
((v[(nl)-1]==-322.0)&&(v[(nh)+1]==-722.0)))
#define valid_dvector_b(v) ( (v!=(double *)NULL) && (*(v)==-322.0) )
#define valid_vector_b(v) ( (v!=(float *)NULL) && (*(v)==-322.0) )
#define valid_matrix(m,nrl,nrh,ncl,nch) ((m!=(float **)NULL) && \
((m[(nrl)-2]==(float*)0x555) \
&&(m[(nrl)][(ncl)-1]==-422.00)&&(m[nrh][nch+1]==-822.0)))
#define valid_dmatrix(m,nrl,nrh,ncl,nch) ((m!=(double **)NULL) && \
((m[(nrl)-2]==(double*)0x555) \
&&(m[(nrl)][(ncl)-1]==-422.00)&&(m[nrh][nch+1]==-822.0)))
#define valid_matrix_b(m) ( (m != (float **)NULL ) && ((*(m-1)==(float*)0x555)&&(*m[1]==-422.00)))
#define valid_dmatrix_b(m) ((m != (double **)NULL ) && ((*(m-1)==(double*)0x555)&&(*m[1]==-422.00)))
#define valid_d4vector( y ) valid_dvector( y, 1, 4 )
#define valid_d4vector_b( y ) valid_dvector_b( y )
#define valid_dframe( y ) valid_dmatrix( y, 1, 4, 1, 4 )
#define NoText (char *)NULL
#define NoFile (FILE *)NULL
/* Global storage */
static linenum = 0;
static char *filetext = NoText;
static char *funcname = NoText;
static error_connect = 0;
static char ms[50];
/* function declaration */
void nrerror_handler( void(*handler)(char error_text[]) ); /* nrlin.c */
/* Error handler derived from Numerical Recipes standard error handler */
void nrerror_DNR(char error_text[])
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
if ( linenum != 0 ){
fprintf(stderr,"Call to %s at line %d in %s\n",
funcname, linenum, filetext);
} else {
fprintf(stderr,"Call location unknown\n");
}
fprintf(stderr,"Press <ENTER> to return to system...");
getchar();
exit(1);}
/* Wrapper functions call this to save information about caller */
void RegisterCall_DNR( int lineref, char *fileref, char *func )
{
if ( !error_connect ) nrerror_handler(nrerror_DNR);
error_connect = 1;
linenum = lineref;
filetext = fileref;
funcname = func;
}
/* Wrapper functions call this to remove information about caller */
void UnRegisterCall_DNR( void )
{
linenum = 0;
}
/* Check out a text string parameter */
void Text_DNR( char *text, int argnum )
{
int i, n;
unsigned char *textp;
if ( text == NoText ) {
sprintf( ms, "NULL text pointer (arg %d)", argnum);
nrerror_DNR( ms );
} else {
n = strlen(text);
if ( n > 255 ) {
sprintf( ms, "text string > 255 chars (arg %d)", argnum);
nrerror_DNR( ms );
} else {
if ( n > 10 ) n = 10;
textp = (unsigned char *)text;
for (i=0; i<n; i++) {
if ( textp[i] < ' ' || textp[i] > 127 ) {
if ( textp[i] != '\n' && textp[i] != '\t' ) {
sprintf( ms, "text string contains strange chars (arg %d)",
argnum);
nrerror_DNR( ms );
}
}
}
}
}
}
/* Check out a 4*4 dmatrix parameter */
void dmat4_DNR( double **a, int argnum )
{
if ( !valid_dmatrix( a, 1, 4, 1, 4 ) ) {
sprintf( ms, "arg %d is not a valid 4*4 dmatrix", argnum);
nrerror_DNR( ms );
}
}
/* Check out a dmatrix parameter */
void dmat_DNR( double **a, int argnum )
{
if ( !valid_dmatrix_b( a ) ) {
sprintf( ms, "arg %d is not a valid dmatrix", argnum);
nrerror_DNR( ms );
}
}
/* Check out a dvector parameter */
void dvec_DNR( double *a, int argnum )
{
if ( !valid_dvector_b( a ) ) { sprintf( ms, "arg %d is not a valid dvector", argnum);
nrerror_DNR( ms );
}
}
void mat_DNR( float **a, int argnum )
{
if ( !valid_matrix_b( a ) ) {
sprintf( ms, "arg %d is not a valid matrix", argnum);
nrerror_DNR( ms );
}
}
void vec_DNR( float *a, int argnum )
{
if ( !valid_vector_b( a ) ) {
sprintf( ms, "arg %d is not a valid vector", argnum);
nrerror_DNR( ms );
}
}
/* Check out a file parameter */
void File_DNR( FILE *fileptr, int argnum )
{
int i;
if ( fileptr == NoFile ) {
sprintf( ms, "NULL file pointer (arg %d)", argnum);
nrerror_DNR( ms );
} else {
i = ferror(fileptr);
if (i) {
sprintf( ms, "file error %d(?) (arg %d)", argnum, argnum);
nrerror_DNR( ms );
}
}
}
/* WRAPPER FUNCTIONS BEGIN HERE */
void free_dmatrix_DNR( int lineref, char *fileref,
double **a, int nrl, int nrh, int ncl, int nch )
{
RegisterCall_DNR( lineref, fileref, "free_dmatrix" );
dmat_DNR( a, 1);
free_dmatrix( a, nrl, nrh, ncl, nch );
UnRegisterCall_DNR();
}
void free_dvector_DNR( int lineref, char *fileref,
double *a, int nrl, int nrh )
{
RegisterCall_DNR( lineref, fileref, "free_dvector" );
dvec_DNR( a, 1);
free_dvector( a, nrl, nrh );
UnRegisterCall_DNR();
}
double **dmatrix_DNR( int lineref, char *fileref,
int nrl, int nrh, int ncl, int nch )
{
double **y;
RegisterCall_DNR( lineref, fileref, "dmatrix" );
y = dmatrix( nrl, nrh, ncl, nch );
UnRegisterCall_DNR();
return(y);
}
double *dvector_DNR( int lineref, char *fileref,
int ncl, int nch )
{
double *y;
RegisterCall_DNR( lineref, fileref, "dvector" );
y = dvector( ncl, nch );
UnRegisterCall_DNR();
return(y);
}
void dmdump_DNR( int lineref, char *fileref,
FILE *outf, char *text, double **a, int a_rows, int a_cols, char *format)
{
RegisterCall_DNR( lineref, fileref, "dmdump" );
File_DNR( outf, 1 );
Text_DNR( text, 2 );
Text_DNR( format, 6);
dmat_DNR( a, 3);
dmdump( outf, text, a, a_rows, a_cols, format);
UnRegisterCall_DNR();
}
void dWriteMatrix_DNR( int lineref, char *fileref,
FILE *outf, char *text, double **a, int a_rows, int a_cols, char *format)
{
RegisterCall_DNR( lineref, fileref, "dWriteMatrix" );
File_DNR( outf, 1 );
Text_DNR( text, 2 );
Text_DNR( format, 6);
dmat_DNR( a, 3);
dWriteMatrix( outf, text, a, a_rows, a_cols, format);
UnRegisterCall_DNR();
}
void dvdump_DNR( int lineref, char *fileref,
FILE *outf, char *text, double *a, int a_els, char *format)
{
RegisterCall_DNR( lineref, fileref, "dvdump" );
File_DNR( outf, 1 );
Text_DNR( text, 2 );
Text_DNR( format, 5);
dvec_DNR( a, 3);
dvdump( outf, text, a, a_els, format);
UnRegisterCall_DNR();
}
void dWriteVector_DNR( int lineref, char *fileref,
FILE *outf, char *text, double *a, int a_els, char *format)
{
RegisterCall_DNR( lineref, fileref, "dWriteVector" );
File_DNR( outf, 1 );
Text_DNR( text, 2 );
Text_DNR( format, 5);
dvec_DNR( a, 3);
dWriteVector( outf, text, a, a_els, format);
UnRegisterCall_DNR();
}
void dWriteScalar_DNR( int lineref, char *fileref,
FILE *outf, char *text, double a, char *format)
{
RegisterCall_DNR( lineref, fileref, "dWriteScalar" );
File_DNR( outf, 1 );
Text_DNR( text, 2 );
Text_DNR( format, 4);
dWriteScalar( outf, text, a, format);
UnRegisterCall_DNR();
}
void dReadMatrix_DNR( int lineref, char *fileref,
FILE *datafile, double **y, int *rows, int *cols, int *error)
/* read a matrix _DNR( int lineref, char *fileref,
NR type) */
{
RegisterCall_DNR( lineref, fileref, "dReadMatrix" );
File_DNR( datafile, 1 );
dmat_DNR( y, 2);
dReadMatrix( datafile, y, rows, cols, error);
UnRegisterCall_DNR();
}
void dReadVector_DNR( int lineref, char *fileref,
FILE *datafile, double *y, int *els, int *error)
/* read a vector (NR type) */
{
RegisterCall_DNR( lineref, fileref, "dReadVector" );
File_DNR( datafile, 1 );
dvec_DNR( y, 2);
dReadVector( datafile, y, els, error);
UnRegisterCall_DNR();
}
double dReadScalar_DNR( int lineref, char *fileref,
FILE *datafile, int *error)
{
double y;
RegisterCall_DNR( lineref, fileref, "dReadScalar" );
File_DNR( datafile, 1 );
y = dReadScalar( datafile, error);
UnRegisterCall_DNR();
return(y);
}
void dmmult_DNR( int lineref, char *fileref,
double **a, int a_rows, int a_cols,
double **b, int b_rows, int b_cols, double **y)
/* multiply two matrices a, b, result in y. y must not be same as a or b */
{
RegisterCall_DNR( lineref, fileref, "dmmult" );
dmat_DNR( a, 1 );
dmat_DNR( b, 4 );
dmat_DNR( y, 7 );
dmmult( a, a_rows, a_cols, b, b_rows, b_cols, y );
UnRegisterCall_DNR();
}
void dmvmult_DNR( int lineref, char *fileref,
double **a, int a_rows, int a_cols, double *b, int b_els, double *y)
/* multiply a matrix a by vector b, result in y. y can be same as b */
{
RegisterCall_DNR( lineref, fileref, "dmvmult" );
dmat_DNR( a, 1 );
dvec_DNR( b, 4 );
dvec_DNR( y, 6 );
dmvmult( a, a_rows, a_cols, b, b_els, y );
UnRegisterCall_DNR();
}
void dmadd_DNR( int lineref, char *fileref,
double **a, int a_rows, int a_cols, double **b, double **y)
/* add two matrices a, b, result in y. y can be same as a or b */
{
RegisterCall_DNR( lineref, fileref, "dmadd" );
dmat_DNR( a, 1 );
dmat_DNR( b, 4 );
dmat_DNR( y, 5 );
dmadd( a, a_rows, a_cols, b, y );
UnRegisterCall_DNR();
}
void dmsmy_DNR( int lineref, char *fileref,
double **a, int a_rows, int a_cols, double r, double **y)
/* multiply a by scalar r, result in y. y can be same as a */
{
RegisterCall_DNR( lineref, fileref, "dmsmy" );
dmat_DNR( a, 1 );
dmat_DNR( y, 5 );
dmsmy( a, a_rows, a_cols, r, y );
UnRegisterCall_DNR();
}
void dmsub_DNR( int lineref, char *fileref,
double **a, int a_rows, int a_cols, double **b, double **y)
/* subtract two matrices a, b, result in y. y can be same as a or b */
{
RegisterCall_DNR( lineref, fileref, "dmsub" );
dmat_DNR( a, 1 );
dmat_DNR( b, 4 );
dmat_DNR( y, 5 );
dmsub( a, a_rows, a_cols, b, y );
UnRegisterCall_DNR();
}
void dmtranspose_DNR( int lineref, char *fileref,
double **a, int a_rows, int a_cols, double **y)
/* transpose matrix a, result in y. y must not be same as a */
{
RegisterCall_DNR( lineref, fileref, "dmtranspose" );
dmat_DNR( a, 1 );
dmat_DNR( y, 4 );
dmtranspose( a, a_rows, a_cols, y );
UnRegisterCall_DNR();
}
void dmfillUT_DNR( int lineref, char *fileref,
double **a, int a_rows, int a_cols)
/* fill upper half of a (square) (lower triangular form) to make a symmetric
if not square, will do the best possible */
{
RegisterCall_DNR( lineref, fileref, "dmfillUT" );
dmat_DNR( a, 1 );
dmfillUT( a, a_rows, a_cols);
UnRegisterCall_DNR();
}
void dvadd_DNR( int lineref, char *fileref,
double *a, int a_els, double *b, double *y)
{
RegisterCall_DNR( lineref, fileref, "dvadd" );
dvec_DNR( a, 1 );
dvec_DNR( b, 3 );
dvec_DNR( y, 4 );
dvadd( a, a_els, b, y);
UnRegisterCall_DNR();
}
void dvsub_DNR( int lineref, char *fileref,
double *a, int a_els, double *b, double *y)
{
RegisterCall_DNR( lineref, fileref, "dvsub" );
dvec_DNR( a, 1 );
dvec_DNR( b, 3 );
dvec_DNR( y, 4 );
dvsub( a, a_els, b, y);
UnRegisterCall_DNR();
}
double dvdot_DNR( int lineref, char *fileref,
double *a, int a_els, double *b){
double y;
RegisterCall_DNR( lineref, fileref, "dvdot" );
dvec_DNR( a, 1 );
dvec_DNR( b, 3 );
y = dvdot( a, a_els, b);
UnRegisterCall_DNR();
return(y);
}
double dvmag_DNR( int lineref, char *fileref,
double *a, int a_els)
{
double y;
RegisterCall_DNR( lineref, fileref, "dvmag" );
dvec_DNR( a, 1 );
y = dvmag( a, a_els);
UnRegisterCall_DNR();
return(y);
}
double dvmnorm_DNR( int lineref, char *fileref,
double *a, int a_els)
{
double y;
RegisterCall_DNR( lineref, fileref, "dvmnorm" );
dvec_DNR( a, 1 );
y = dvmnorm( a, a_els);
UnRegisterCall_DNR();
return(y);
}
void dvsmy_DNR( int lineref, char *fileref,
double *a, int a_els, double r, double *y)
{
RegisterCall_DNR( lineref, fileref, "dvsmy" );
dvec_DNR( a, 1 );
dvec_DNR( y, 4 );
dvsmy( a, a_els, r, y);
UnRegisterCall_DNR();
}
void dvpiv_DNR( int lineref, char *fileref,
double *a, int a_els, double r, double *b, double *y)
{
RegisterCall_DNR( lineref, fileref, "dvpiv" );
dvec_DNR( a, 1 );
dvec_DNR( b, 4 );
dvec_DNR( y, 5 );
dvpiv( a, a_els, r, b, y);
UnRegisterCall_DNR();
}
void dmcopy_DNR( int lineref, char *fileref,
double **a, int a_rows, int a_cols, double **b) /* copy a matrix */
{
RegisterCall_DNR( lineref, fileref, "dmcopy" );
dmat_DNR( a, 1 );
dmat_DNR( b, 4 );
dmcopy( a, a_rows, a_cols, b);
UnRegisterCall_DNR();
}
void dvcopy_DNR( int lineref, char *fileref,
double *a, int a_els, double *y) /* copy a vector */{
RegisterCall_DNR( lineref, fileref, "dvcopy" );
dvec_DNR( a, 1 );
dvec_DNR( y, 3 );
dvcopy( a, a_els, y);
UnRegisterCall_DNR();
}
double dvcomp_DNR( int lineref, char *fileref,
double *a, int a_els, double *b) /* compare two vectors */
{
double y;
RegisterCall_DNR( lineref, fileref, "dvcomp" );
dvec_DNR( a, 1 );
dvec_DNR( b, 3 );
y = dvcomp( a, a_els, b);
UnRegisterCall_DNR();
return(y);
}
void dgetcolumn_DNR( int lineref, char *fileref,
double **G, int col, double *v, int nels) /* get a column from a matrix */
{
RegisterCall_DNR( lineref, fileref, "dgetcolumn" );
dmat_DNR( G, 1 );
dvec_DNR( v, 3 );
dgetcolumn( G, col, v, nels );
UnRegisterCall_DNR();
}
void dputcolumn_DNR( int lineref, char *fileref,
double *v, int nels, double **G, int col) /* put a column into a matrix */
{
RegisterCall_DNR( lineref, fileref, "dputcolumn" );
dmat_DNR( G, 3 );
dvec_DNR( v, 1 );
dputcolumn( v, nels, G, col );
UnRegisterCall_DNR();
}
void dinverse_DNR( int lineref, char *fileref,
double **a, int n, double **y )
{
RegisterCall_DNR( lineref, fileref, "dinverse" );
dmat_DNR( a, 1 );
dmat_DNR( y, 3 );
dinverse( a, n, y );
UnRegisterCall_DNR();
}
void dinverse_mult_DNR( int lineref, char *fileref,
double **a, int a_rows, double **b, int b_cols, double **y )
{
RegisterCall_DNR( lineref, fileref, "dinverse_mult" );
dmat_DNR( a, 1 );
dmat_DNR( b, 3 );
dmat_DNR( y, 5 );
dinverse_mult( a, a_rows, b, b_cols, y);
UnRegisterCall_DNR();
}
void dPDSinverse_DNR( int lineref, char *fileref,
double **a, int n, double **y )
{
RegisterCall_DNR( lineref, fileref, "dPDSinverse" );
dmat_DNR( a, 1 );
dmat_DNR( y, 3 );
dPDSinverse(a, n, y ); UnRegisterCall_DNR();
}
void dPDS_L_inverse_DNR( int lineref, char *fileref,
double **a, int n, double **y )
{
RegisterCall_DNR( lineref, fileref, "dPDS_L_inverse" );
dmat_DNR( a, 1 );
dmat_DNR( y, 3 );
dPDS_L_inverse( a, n, y );
UnRegisterCall_DNR();
}
void choldc_DNR( int lineref, char *fileref,
float **a, int n, float p[])
{
RegisterCall_DNR( lineref, fileref, "choldc" );
mat_DNR( a, 1 );
vec_DNR( p, 3 );
choldc( a, n, p );
UnRegisterCall_DNR();
}
void cholsl_DNR( int lineref, char *fileref,
float **a, int n, float p[], float b[], float x[])
{
RegisterCall_DNR( lineref, fileref, "cholsl" );
mat_DNR( a, 1 );
vec_DNR( p, 3 );
vec_DNR( b, 4 );
vec_DNR( x, 5 );
cholsl( a, n, p, b, x );
UnRegisterCall_DNR();
}
void lubksb_DNR( int lineref, char *fileref,
float **a, int n, int *indx, float b[])
{
RegisterCall_DNR( lineref, fileref, "lubksb" );
mat_DNR( a, 1 );
vec_DNR( b, 4 );
lubksb( a, n, indx, b );
UnRegisterCall_DNR();
}
void ludcmp_DNR( int lineref, char *fileref,
float **a, int n, int *indx, float *d)
{
RegisterCall_DNR( lineref, fileref, "ludcmp" );
mat_DNR( a, 1 );
vec_DNR( d, 4 );
ludcmp( a, n, indx, d );
UnRegisterCall_DNR();
}
void dcholdc_DNR( int lineref, char *fileref,
double **a, int n, double p[])
{
RegisterCall_DNR( lineref, fileref, "dcholdc" );
dmat_DNR( a, 3 );
dvec_DNR( p, 3 );
dcholdc( a, n, p );
UnRegisterCall_DNR();}
void dcholsl_DNR( int lineref, char *fileref,
double **a, int n, double p[], double b[], double x[])
{
RegisterCall_DNR( lineref, fileref, "dcholsl" );
dmat_DNR( a, 1 );
dvec_DNR( p, 3 );
dvec_DNR( b, 4 );
dvec_DNR( x, 5 );
dcholsl( a, n, p, b, x );
UnRegisterCall_DNR();
}
void dlubksb_DNR( int lineref, char *fileref,
double **a, int n, int *indx, double b[])
{
RegisterCall_DNR( lineref, fileref, "dlubksb" );
dmat_DNR( a, 1 );
dvec_DNR( b, 4 );
dlubksb( a, n, indx, b );
UnRegisterCall_DNR();
}
void dludcmp_DNR( int lineref, char *fileref,
double **a, int n, int *indx, double *d)
{
RegisterCall_DNR( lineref, fileref, "dludcmp" );
dmat_DNR( a, 1 );
dvec_DNR( d, 4 );
dludcmp( a, n, indx, d );
UnRegisterCall_DNR();
}
double **dframe_DNR( int lineref, char *fileref )
{
double **y;
RegisterCall_DNR( lineref, fileref, "dframe" );
y = dframe();
UnRegisterCall_DNR();
return y;
}
double *d4vector_DNR( int lineref, char *fileref )
{
double *y;
RegisterCall_DNR( lineref, fileref, "d4vector" );
y = d4vector();
UnRegisterCall_DNR();
return y;
}
void free_d4vector_DNR( int lineref, char *fileref, double *y )
{
RegisterCall_DNR( lineref, fileref, "free_d4vector" );
free_d4vector(y);
UnRegisterCall_DNR();
}
void free_dframe_DNR( int lineref, char *fileref, double **y )
{ RegisterCall_DNR( lineref, fileref, "free_d4vector" );
free_dframe(y);
UnRegisterCall_DNR();
}
void chainmult_DNR( int lineref, char *fileref,
double **a, double **b )
{
RegisterCall_DNR( lineref, fileref, "chainmult" );
dmat4_DNR( a, 1 );
dmat4_DNR( b, 2 );
chainmult( a, b );
UnRegisterCall_DNR();
}
void rotframe_DNR( int lineref, char *fileref,
double **a, int axis, double theta)
{
RegisterCall_DNR( lineref, fileref, "rotframe" );
dmat4_DNR( a, 1 );
rotframe( a, axis, theta );
UnRegisterCall_DNR();
}
void transframe_DNR( int lineref, char *fileref,
double **a, double dx, double dy, double dz)
{
RegisterCall_DNR( lineref, fileref, "transframe" );
dmat4_DNR( a, 1 );
transframe( a, dx, dy, dz );
UnRegisterCall_DNR();
}
void orthog_DNR( int lineref, char *fileref,
double **a )
{
RegisterCall_DNR( lineref, fileref, "orthog" );
dmat4_DNR( a, 1 );
orthog( a );
UnRegisterCall_DNR();
}
void orthogKI_DNR( int lineref, char *fileref,
double **a )
{
RegisterCall_DNR( lineref, fileref, "orthogKI" );
dmat4_DNR( a, 1 );
orthogKI( a );
UnRegisterCall_DNR();
}
void dframeinverse_DNR( int lineref, char *fileref,
double **E, double **Ein )
{
RegisterCall_DNR( lineref, fileref, "dframeinverse" );
dmat4_DNR( E, 1 );
dmat4_DNR( Ein, 2 );
dframeinverse( E, Ein );
UnRegisterCall_DNR();
}
void dcross4_DNR( int lineref, char *fileref,
double *a, double *b, double *c ){
RegisterCall_DNR( lineref, fileref, "dcross4" );
dvec_DNR( a, 1 );
dvec_DNR( b, 2 );
dvec_DNR( c, 3 );
dcross4( a, b, c );
UnRegisterCall_DNR();
}
double column_dot_DNR( int lineref, char *fileref,
double **G1, int col1, double **G2, int col2 )
{
double y;
RegisterCall_DNR( lineref, fileref, "dcolumn_dot" );
dmat4_DNR( G1, 1 );
dmat4_DNR( G2, 3 );
y = column_dot( G1, col1, G2, col2 );
UnRegisterCall_DNR();
return (y);
}
void make_d_omega_DNR( int lineref, char *fileref,
double **G1, double **G2, double *w)
{
RegisterCall_DNR( lineref, fileref, "make_d_omega" );
dmat4_DNR( G1, 1 );
dmat4_DNR( G2, 2 );
make_d_omega( G1, G2, w );
UnRegisterCall_DNR();
}
void compare_link_frames_DNR( int lineref, char *fileref,
double **G1, double **G2, double *wp, double *wr, double *ww )
{
RegisterCall_DNR( lineref, fileref, "compare_link_frames" );
dmat4_DNR( G1, 1 );
dmat4_DNR( G2, 2 );
compare_link_frames( G1, G2, wp, wr, ww );
UnRegisterCall_DNR();
}
void compare_frames_DNR( int lineref, char *fileref,
double **G1, double **G2, double *wp, double *wr, double *ww )
{
RegisterCall_DNR( lineref, fileref, "compare_frames" );
dmat4_DNR( G1, 1 );
dmat4_DNR( G2, 2 );
compare_frames( G1, G2, wp, wr, ww );
UnRegisterCall_DNR();
}
void FindAxis_DNR( int lineref, char *fileref,
double **T1, double **T2, double *axis,
double *sphi, double *cphi,int *twitching)
{
RegisterCall_DNR( lineref, fileref, "FindAxis" );
dmat4_DNR( T1, 1 );
dmat4_DNR( T2, 2 );
dvec_DNR( axis, 3 );
FindAxis( T1, T2, axis, sphi, cphi, twitching );
UnRegisterCall_DNR();
}
void RotateVector_DNR( int lineref, char *fileref,
double *v, double *axis, double sphi, double cphi, double *y)
{
RegisterCall_DNR( lineref, fileref, "RotateVector" );
dvec_DNR( v, 1 );
dvec_DNR( axis, 2 );
dvec_DNR( y, 5 );
RotateVector( v, axis, sphi, cphi, y );
UnRegisterCall_DNR();
}
void RotateFrame_DNR( int lineref, char *fileref,
double **T1, double *axis, double sphi, double cphi, double **T2)
{
RegisterCall_DNR( lineref, fileref, "RotateFrame" );
dmat_DNR( T1, 1 );
dvec_DNR( axis, 2 );
dmat_DNR( T2, 5 );
RotateFrame( T1, axis, sphi, cphi, T2 );
UnRegisterCall_DNR();
}
#ifdef TEST
#include "nrcheck.h"
void main( void )
{
/* test program for above utility routines */
double **a, **b, **c, **bT;
double *x, *y, *z;
FILE *infile, *outfile;
int a_rows, a_cols, b_rows, b_cols, errors, xn, yn;
char L1[] = {"Matrix A"};
char form[] = {"%11.4lf"};
infile = fopen("mat.dat", "r");
outfile = fopen("mat.out", "w");
a = dmatrix(1,10,1,10);
b = dmatrix(1,10,1,10);
x = dvector(1,10);
y = dvector(1,10);
dReadMatrix( infile, a, &a_rows, &a_cols, &errors);
printf("Reading a done - hit <return> to go on\n");
getchar();
dReadMatrix( infile, b, &b_rows, &b_cols, &errors);
dReadVector( infile, x, &xn, &errors);
dReadVector( infile, y, &yn, &errors);
printf("Reading done - hit <return> to go on\n");
getchar();
dmdump( stderr, L1, a, a_rows, a_cols, form);
dmdump( outfile, "Matrix B", b, b_rows, b_cols, "%8.2lf");
dvdump( outfile, "Vector x", x, xn, form);
dvdump( outfile, "Vector y", y, yn, "%8.2lf");
z = dvector( 1, xn );
dvadd( x, xn, y, z );
dvdump( outfile, "x + y", z, xn, "%8.2lf");
dvsub( x, xn, y, z );
dvdump( outfile, "x - y", z, xn, "%8.2lf");
dvsmy( x, xn, 2.0, z );
dvdump( outfile, "2x", z, xn, "%8.2lf");
printf("Magnitude of 2x: %7.2lf\n", dvmag( z, xn ));
printf("dot product x.y: %7.2lf\n", dvdot( x, xn, y));
dmvmult( a, a_rows, a_cols, x, xn, z );
dvdump( outfile, "Ax", z, xn, "%8.2lf"); if ( !valid_dmatrix( a, 1, 10, 1, 10) )nrerror("a invalid");
c = dmatrix( 1, a_rows, 1, b_cols );
bT = dmatrix( 1, b_cols, 1, b_rows );
dmtranspose( b, b_rows, b_cols, bT);
dmdump( outfile, "Matrix B (transposed)", bT, b_cols, b_rows, "%8.2lf");
dmmult( a, a_rows, a_cols, bT, b_cols, b_rows, c);
dmdump( outfile, "Matrix AB", c, a_rows, b_rows, "%8.2lf");
dmdump( outfile, L1, a, a_rows, a_cols, form);
/* dmfillUT( a, a_rows, a_cols );
dmdump( outfile, "Symmetrified matrix A", a, a_rows, a_cols, "%8.2lf"); */
reportmemory(stdout);
free_dmatrix( a, 1, 10, 1, 10);
free_dmatrix( b, 1, 10, 1, 10);
free_dmatrix( c, 1, a_rows, 1, b_cols);
free_dmatrix( bT, 1, b_cols, 1, b_rows);
free_dvector( x, 1, 10 );
free_dvector( y, 1, 10 );
free_dvector( z, 1, xn );
reportmemory(stdout);
printf("Type <return> to exit\n");
getchar();
}
#endif