/* //: var1.c -- vector AR(1) data generator */ /* Program Usage and Configuration: */ const char *pstUsage="\ var1 [-etvh?] [-n yyyy] [-k zzzz] [-s wwww] [-i ]\n\ \n\ This program generates multivariate AR(1) based on the settings in \n\ the configuration file, var1.ini by default. \n\ If the EDITOR environmental variable is set, the -e option will edit\n\ a the configuration file, or create a sample one if it does not exist.\n\ Options:\n\ -e Creates or edits the configuration file using environment variable 'EDITOR'\n\ -h or -? this information\n\ -n number of total observations\n\ -k number of in control observations\n\ -s random number seed\n\ -i uses this as a configuration file \n\ -v verbose: print out the setup information\n\ -t terse: #\n\ \n\ Dave Forrest 7 July 1999\n\ "; const char *pstIniDefault ="\ # var1 configuration file\n\ # \n\ # Multivariate VAR1 initialization file\n\ # these are the configurable parameters:\n\ # this file set up a multivariate AR(1) process\n\ # with the following parameters\n\ \n\ # Process Mean\n\ mean = 0 0 # mean data vector\n\ \n\ # Innovation or shift\n\ shift = 1 0 # shift in the uncorrelated noise process\n\ \n\ # covariance of the noise process\n\ sigma = \"1.00 1.10\n\ 1.10 2.00\" # Covariance of data vectors\n\ \n\ # Autoregressive Parameter\n\ phi = \"0.50 0.00\n\ 0.00 -0.25\"\n\ \n\ # These parameters may be overridden from the command line \n\ # total number of observations (option -n xxxx) \n\ n_obs = 15 \n\ # number of in control points (option -k yyyy)\n\ n_incontrol = 3 \n\ # random number seed (option -s zzzz)\n\ seed = 0 \n\ \n\ #Also, try var1 -? for information on the other options:\n\ # -v verbose\n\ #\n\ # Source code developed on drf5n@hob.private:src/var1_1/var1.c\n\ # David Forrest drf5n@virginia.edu 2000-01-18\n\ # http://mug.sys.virginia.edu/~drf5n/\n\ "; /* compile with: Simple Compilation, which seems to work on linux, AIX, DEC: cc var1.c -lm -o var1 For the Sun machine which uses ucbcc, append -Xc to eliminate the unix extentions: cc var1.c -lm -o var1 -Xc I used gnu's gcc with the following options during development: gcc -Wall -pedantic var1.c -lm -o var1 Microsoft Visual C: Use the menus to build var1.exe */ /* in Microsoft environments, the unistd.h file and its attendant functions do not exist, so, add a define for Microsoft environments */ #ifndef ANSI_C #define ANSI_C #endif #ifndef DEFAULT_EDITOR # if defined unix | defined __unix__ | defined __unix # define DEFAULT_EDITOR "vi" # else # define DEFAULT_EDITOR "notepad" # endif #endif #ifndef PROTOTYPES_IN_STRUCT #define PROTOTYPES_IN_STRUCT #endif #define DEBUG 0 #if DEBUG > 5 #define DebugCode( code ) { code;}; #else #define DebugCode( code ) /**/ #endif #include #include #include #include #include /* Declarations for getopt. Copyright (C) 1989,90,91,92,93,94,96,97 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. The GNU C Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more details. You should have received a copy of the GNU Library General Public License along with the GNU C Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #ifndef _GETOPT_H #define _GETOPT_H 1 #ifdef __cplusplus extern "C" { #endif /* For communication from `getopt' to the caller. When `getopt' finds an option that takes an argument, the argument value is returned here. Also, when `ordering' is RETURN_IN_ORDER, each non-option ARGV-element is returned here. */ char *optarg; /* Index in ARGV of the next element to be scanned. This is used for communication to and from the caller and for communication between successive calls to `getopt'. On entry to `getopt', zero means this is the first call; initialize. When `getopt' returns -1, this is the index of the first of the non-option elements that the caller should itself scan. Otherwise, `optind' communicates from one call to the next how much of ARGV has been scanned so far. */ extern int optind; /* Callers store zero here to inhibit the error message `getopt' prints for unrecognized options. */ extern int opterr; /* Set to an option character which was unrecognized. */ extern int optopt; /* Describe the long-named options requested by the application. The LONG_OPTIONS argument to getopt_long or getopt_long_only is a vector of `struct option' terminated by an element containing a name which is zero. The field `has_arg' is: no_argument (or 0) if the option does not take an argument, required_argument (or 1) if the option requires an argument, optional_argument (or 2) if the option takes an optional argument. If the field `flag' is not NULL, it points to a variable that is set to the value given in the field `val' when the option is found, but left unchanged if the option is not found. To have a long-named option do something other than set an `int' to a compiled-in constant, such as set a value from `optarg', set the option's `flag' field to zero and its `val' field to a nonzero value (the equivalent single-letter option character, if there is one). For long options that have a zero `flag' field, `getopt' returns the contents of the `val' field. */ struct option { #if defined (__STDC__) && __STDC__ const char *name; #else char *name; #endif /* has_arg can't be an enum because some compilers complain about type mismatches in all the code that assumes it is an int. */ int has_arg; int *flag; int val; }; /* Names for the values of the `has_arg' field of `struct option'. */ #define no_argument 0 #define required_argument 1 #define optional_argument 2 #if defined (__STDC__) && __STDC__ #ifdef __GNU_LIBRARY__ /* Many other libraries have conflicting prototypes for getopt, with differences in the consts, in stdlib.h. To avoid compilation errors, only prototype getopt for the GNU C library. */ extern int getopt (int argc, char *const *argv, const char *shortopts); #else /* not __GNU_LIBRARY__ */ extern int getopt (); #endif /* __GNU_LIBRARY__ */ extern int getopt_long (int argc, char *const *argv, const char *shortopts, const struct option *longopts, int *longind); extern int getopt_long_only (int argc, char *const *argv, const char *shortopts, const struct option *longopts, int *longind); /* Internal only. Users should not call this directly. */ extern int _getopt_internal (int argc, char *const *argv, const char *shortopts, const struct option *longopts, int *longind, int long_only); #else /* not __STDC__ */ extern int getopt (); extern int getopt_long (); extern int getopt_long_only (); extern int _getopt_internal (); #endif /* __STDC__ */ #ifdef __cplusplus } #endif #endif /* _GETOPT_H */ #define isatty(X) (0) /* make the matrix libraries not do terminal IO */ /* included header files for routines */ /************************************************************************** ** ** Copyright (C) 1993 David E. Stewart & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* err.h 28/09/1993 */ /* RCS id: $Id: var1.c,v 1.20 2000/02/01 21:11:27 drf5n Exp drf5n $ */ #ifndef ERRHEADER #define ERRHEADER #include /* #include "machine.h" */ /* Error recovery */ extern jmp_buf restart; /* max. # of error lists */ #define ERR_LIST_MAX_LEN 10 /* main error functions */ #ifndef ANSI_C extern int ev_err(); /* main error handler */ extern int set_err_flag(); /* for different ways of handling errors, returns old value */ extern int count_errs(); /* to avoid "too many errors" */ extern int err_list_attach(); /* for attaching a list of errors */ extern int err_is_list_attached(); /* checking if a list is attached */ extern int err_list_free(); /* freeing a list of errors */ #else /* ANSI_C */ extern int ev_err(char *,int,int,char *,int); /* main error handler */ extern int set_err_flag(int flag); /* for different ways of handling errors, returns old value */ extern int count_errs(int true_false); /* to avoid "too many errors" */ extern int err_list_attach(int list_num, int list_len, char **err_ptr,int warn); /* for attaching a list of errors */ extern int err_is_list_attached(int list_num); /* checking if a list is attached */ extern int err_list_free(int list_num); /* freeing a list of errors */ #endif /* error(E_TYPE,"myfunc") raises error type E_TYPE for function my_func() */ #define error(err_num,fn_name) ev_err(__FILE__,err_num,__LINE__,fn_name,0) /* warning(WARN_TYPE,"myfunc") raises warning type WARN_TYPE for function my_func() */ #define warning(err_num,fn_name) ev_err(__FILE__,err_num,__LINE__,fn_name,1) /* error flags */ #define EF_EXIT 0 /* exit on error */ #define EF_ABORT 1 /* abort (dump core) on error */ #define EF_JUMP 2 /* jump on error */ #define EF_SILENT 3 /* jump, but don't print message */ #define ERREXIT() set_err_flag(EF_EXIT) #define ERRABORT() set_err_flag(EF_ABORT) /* don't print message */ #define SILENTERR() if ( ! setjmp(restart) ) set_err_flag(EF_SILENT) /* return here on error */ #define ON_ERROR() if ( ! setjmp(restart) ) set_err_flag(EF_JUMP) /* error types */ #define E_UNKNOWN 0 #define E_SIZES 1 #define E_BOUNDS 2 #define E_MEM 3 #define E_SING 4 #define E_POSDEF 5 #define E_FORMAT 6 #define E_INPUT 7 #define E_NULL 8 #define E_SQUARE 9 #define E_RANGE 10 #define E_INSITU2 11 #define E_INSITU 12 #define E_ITER 13 #define E_CONV 14 #define E_START 15 #define E_SIGNAL 16 #define E_INTERN 17 #define E_EOF 18 #define E_SHARED_VECS 19 #define E_NEG 20 #define E_OVERWRITE 21 #define E_BREAKDOWN 22 /* warning types */ #define WARN_UNKNOWN 0 #define WARN_WRONG_TYPE 1 #define WARN_NO_MARK 2 #define WARN_RES_LESS_0 3 #define WARN_SHARED_VEC 4 /* error catching macros */ /* execute err_part if error errnum is raised while executing ok_part */ #define catch(errnum,ok_part,err_part) \ { jmp_buf _save; int _err_num, _old_flag; \ _old_flag = set_err_flag(EF_SILENT); \ MEM_COPY(restart,_save,sizeof(jmp_buf)); \ if ( (_err_num=setjmp(restart)) == 0 ) \ { ok_part; \ set_err_flag(_old_flag); \ MEM_COPY(_save,restart,sizeof(jmp_buf)); } \ else if ( _err_num == errnum ) \ { set_err_flag(_old_flag); \ MEM_COPY(_save,restart,sizeof(jmp_buf)); \ err_part; } \ else { set_err_flag(_old_flag); \ MEM_COPY(_save,restart,sizeof(jmp_buf)); \ error(_err_num,"catch"); \ } \ } /* execute err_part if any error raised while executing ok_part */ #define catchall(ok_part,err_part) \ { jmp_buf _save; int _err_num, _old_flag; \ _old_flag = set_err_flag(EF_SILENT); \ MEM_COPY(restart,_save,sizeof(jmp_buf)); \ if ( (_err_num=setjmp(restart)) == 0 ) \ { ok_part; \ set_err_flag(_old_flag); \ MEM_COPY(_save,restart,sizeof(jmp_buf)); } \ else \ { set_err_flag(_old_flag); \ MEM_COPY(_save,restart,sizeof(jmp_buf)); \ err_part; } \ } /* print message if error raised while executing ok_part, then re-raise error to trace calls */ #define tracecatch(ok_part,function) \ { jmp_buf _save; int _err_num, _old_flag; \ _old_flag = set_err_flag(EF_JUMP); \ MEM_COPY(restart,_save,sizeof(jmp_buf)); \ if ( (_err_num=setjmp(restart)) == 0 ) \ { ok_part; \ set_err_flag(_old_flag); \ MEM_COPY(_save,restart,sizeof(jmp_buf)); } \ else \ { set_err_flag(_old_flag); \ MEM_COPY(_save,restart,sizeof(jmp_buf)); \ error(_err_num,function); } \ } #endif /* ERRHEADER */ /* * readfloat.h -- Vector reading routines * */ #ifdef _READFLOATS_H_ #else #define _READFLOATS_H_ #include #include #include int readfloats(char *pst, float ** ppvcX, int * pscDimMax); #endif /* * getdelim.h -- this function copies the GNU glibc function * int getdelim( char **ppline, size_t* pscLen, int sc, FILE* pF); * which reads a malloc limited string into the buffer * * Pre condition: * *ppline == NULL && *psclen == 0 : call malloc for new memory * *ppline != NULL && *psclen >= 0 : realloc memory * sc = '/n' or some other delimiter character * pf = stdin or some other file pointer * * Post condition: * Returns number of characters read or EOF * if it cannot find a delimiter, it reads until EOF, and returns the * number of characters read, and an EOF on the next read. * *pline = pointer to null delimited string / read buffer * *psclen = new length of read buffer memory * * dave forrest 1998 July 13 */ int getdelim( char **ppline, size_t* pscLen, int sc, FILE* pF); /* +++Date last modified: 05-Jul-1997 */ /* ======================================================================= CFG.h Configuration file handler. A. Reitsma, Delft, The Netherlands. v1.00 94-07-09 Public Domain. V1.10 97-08-01 drf5n@virginia.edu -- added defaults ---------------------------------------------------------------------- */ struct CfgStrings { char * name ; /* LHS name */ char * data ; /* RHS string */ char * def ; /* Default value */ }; int CfgRead( char * Filename, struct CfgStrings * CfgInfo ); /* ==== CFG.h end ===================================================== */ /*****************************************************/ #define BUFFER_SIZE 1024*16 #define FALSE 0 #define TRUE 1 /* drf5n: POSIX function for terminal/batch identification */ /* extern int fileno __P((FILE*)); */ #ifndef isascii #define isascii(i) ((i&0x7f)==(i)) #endif /***************************************************** include files for the meschach matrix library functions */ /* Configuration for Linux: */ #define U_INT_DEF /* floating point precision */ /* you can choose single, double or long double (if available) precision */ #define FLOAT 1 #define DOUBLE 2 #define LONG_DOUBLE 3 /* #undef REAL_FLT */ /* #undef REAL_DBL */ /* if nothing is defined, choose double precision */ #ifndef REAL_DBL #ifndef REAL_FLT #define REAL_DBL 1 #endif #endif /* single precision */ #ifdef REAL_FLT #define Real float #define LongReal float #define REAL FLOAT #define LONGREAL FLOAT #endif /* double precision */ #ifdef REAL_DBL #define Real double #define LongReal double #define REAL DOUBLE #define LONGREAL DOUBLE #endif /* machine epsilon or unit roundoff error */ /* This is correct on most IEEE Real precision systems */ #ifdef DBL_EPSILON #if REAL == DOUBLE #define MACHEPS DBL_EPSILON #elif REAL == FLOAT #define MACHEPS FLT_EPSILON #elif REAL == LONGDOUBLE #define MACHEPS LDBL_EPSILON #endif #endif #define F_MACHEPS 1.19209e-07 #define D_MACHEPS 2.22045e-16 #ifndef MACHEPS #if REAL == DOUBLE #define MACHEPS D_MACHEPS #elif REAL == FLOAT #define MACHEPS F_MACHEPS #elif REAL == LONGDOUBLE #define MACHEPS D_MACHEPS #endif #endif /* #undef M_MACHEPS */ /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* Type definitions for general purpose maths package */ #ifndef MATRIXH /* RCS id: $Id: var1.c,v 1.20 2000/02/01 21:11:27 drf5n Exp drf5n $ */ #define MATRIXH /* vector definition */ typedef struct { int dim, max_dim; Real *ve; } VEC; /* matrix definition */ typedef struct { int m, n; int max_m, max_n, max_size; Real **me,*base; /* base is base of alloc'd mem */ } MAT; /* band matrix definition */ typedef struct { MAT *mat; /* matrix */ int lb,ub; /* lower and upper bandwidth */ } BAND; /* permutation definition */ typedef struct { int size, max_size, *pe; } PERM; /* integer vector definition */ typedef struct { int dim, max_dim; int *ive; } IVEC; #ifndef MALLOCDECL #ifndef ANSI_C extern char *malloc(), *calloc(), *realloc(); #else extern void *malloc(size_t), *calloc(size_t,size_t), *realloc(void *,size_t); #endif #endif #ifndef ANSI_C extern void m_version(); #else void m_version( void ); #endif #ifndef ANSI_C /* allocate one object of given type */ #define NEW(type) ((type *)calloc(1,sizeof(type))) /* allocate num objects of given type */ #define NEW_A(num,type) ((type *)calloc((unsigned)(num),sizeof(type))) /* re-allocate arry to have num objects of the given type */ #define RENEW(var,num,type) \ ((var)=(type *)((var) ? \ realloc((char *)(var),(unsigned)(num)*sizeof(type)) : \ calloc((unsigned)(num),sizeof(type)))) #define MEMCOPY(from,to,n_items,type) \ MEM_COPY((char *)(from),(char *)(to),(unsigned)(n_items)*sizeof(type)) #else /* allocate one object of given type */ #define NEW(type) ((type *)calloc((size_t)1,(size_t)sizeof(type))) /* allocate num objects of given type */ #define NEW_A(num,type) ((type *)calloc((size_t)(num),(size_t)sizeof(type))) /* re-allocate arry to have num objects of the given type */ #define RENEW(var,num,type) \ ((var)=(type *)((var) ? \ realloc((char *)(var),(size_t)((num)*sizeof(type))) : \ calloc((size_t)(num),(size_t)sizeof(type)))) #define MEMCOPY(from,to,n_items,type) \ MEM_COPY((char *)(from),(char *)(to),(unsigned)(n_items)*sizeof(type)) #endif /* standard copy & zero functions */ #define MEM_COPY(from,to,size) memmove((to),(from),(size)) #define MEM_ZERO(where,size) memset((where),'\0',(size)) /* type independent min and max operations */ #ifndef max #define max(a,b) ((a) > (b) ? (a) : (b)) #endif #ifndef min #define min(a,b) ((a) > (b) ? (b) : (a)) #endif #undef TRUE #define TRUE 1 #undef FALSE #define FALSE 0 /* for input routines */ #define MAXLINE 81 /* Dynamic memory allocation */ /* Should use M_FREE/V_FREE/PX_FREE in programs instead of m/v/px_free() as this is considerably safer -- also provides a simple type check ! */ #ifndef ANSI_C extern VEC *v_get(), *v_resize(); extern MAT *m_get(), *m_resize(); extern PERM *px_get(), *px_resize(); extern IVEC *iv_get(), *iv_resize(); extern int m_free(),v_free(); extern int px_free(); extern int iv_free(); extern BAND *bd_get(), *bd_resize(); extern int bd_free(); #else /* get/resize vector to given dimension */ extern VEC *v_get(int), *v_resize(VEC *,int); /* get/resize matrix to be m x n */ extern MAT *m_get(int,int), *m_resize(MAT *,int,int); /* get/resize permutation to have the given size */ extern PERM *px_get(int), *px_resize(PERM *,int); /* get/resize an integer vector to given dimension */ extern IVEC *iv_get(int), *iv_resize(IVEC *,int); /* get/resize a band matrix to given dimension */ extern BAND *bd_get(int,int,int), *bd_resize(BAND *,int,int,int); /* free (de-allocate) (band) matrices, vectors, permutations and integer vectors */ extern int iv_free(IVEC *); extern int m_free(MAT *),v_free(VEC *),px_free(PERM *); extern int bd_free(BAND *); #endif /* MACROS */ /* macros that also check types and sets pointers to NULL */ #define M_FREE(mat) ( m_free(mat), (mat)=(MAT *)NULL ) #define V_FREE(vec) ( v_free(vec), (vec)=(VEC *)NULL ) #define PX_FREE(px) ( px_free(px), (px)=(PERM *)NULL ) #define IV_FREE(iv) ( iv_free(iv), (iv)=(IVEC *)NULL ) #define MAXDIM 2001 /**************** from meminfo.h */ /* for hash table in mem_stat.c */ /* Note: the hash size should be a prime, or at very least odd */ #define MEM_HASHSIZE 509 #define MEM_HASHSIZE_FILE "meminfo.h" /* default: memory information is off */ /* set it to 1 if you want it all the time */ #define MEM_SWITCH_ON_DEF 0 /* available standard types */ #define TYPE_NULL (-1) #define TYPE_MAT 0 #define TYPE_BAND 1 #define TYPE_PERM 2 #define TYPE_VEC 3 #define TYPE_IVEC 4 #ifdef SPARSE #define TYPE_ITER 5 #define TYPE_SPROW 6 #define TYPE_SPMAT 7 #endif #ifdef COMPLEX #ifdef SPARSE #define TYPE_ZVEC 8 #define TYPE_ZMAT 9 #else #define TYPE_ZVEC 5 #define TYPE_ZMAT 6 #endif #endif /* structure for memory information */ typedef struct { long bytes; /* # of allocated bytes for each type (summary) */ int numvar; /* # of allocated variables for each type */ } MEM_ARRAY; #define mem_info() mem_info_file(stdout,0) #define mem_stat_reg(var,type) mem_stat_reg_list((void **)var,type,0) #define MEM_STAT_REG(var,type) mem_stat_reg_list((void **)&(var),type,0) #define mem_stat_free(mark) mem_stat_free_list(mark,0) #define mem_bytes(type,old_size,new_size) \ mem_bytes_list(type,old_size,new_size,0) #define mem_numvar(type,num) mem_numvar_list(type,num,0) #define TYPE_NULL (-1) #define TYPE_MAT 0 #define TYPE_BAND 1 #define TYPE_PERM 2 #define TYPE_VEC 3 #define TYPE_IVEC 4 /* internal type */ typedef struct { char **type_names; /* array of names of types (strings) */ int (**free_funcs)(void*); /* array of functions for releasing types */ /* array of functions for releasing types cdecl> explain int (*free_funcs[5])(void*) declare free_funcs as array 5 of pointer to function (pointer to void) returning int */ /* int (*free_funcs[5])(void*); */ /* int (**free_funcs)(void*); */ unsigned ntypes; /* max number of types */ MEM_ARRAY *info_sum; /* local array for keeping track of memory */ } MEM_CONNECT; /* max number of lists of types */ #define MEM_CONNECT_MAX_LISTS 5 int mem_stat_reg_list(void **var,int type,int list); int mem_stat_mark(int mark); int mem_stat_free_list(int mark,int list); int mem_stat_show_mark(void); void mem_stat_dump(FILE *fp,int list); /* int mem_attach_list(int list,int ntypes,char *type_names[], int (*free_funcs[])(), MEM_ARRAY info_sum[]); */ int mem_attach_list(int list,int ntypes,char *type_names[], int (*free_funcs[])(void*), MEM_ARRAY info_sum[]); /* Entry level access to data structures */ #ifdef DEBUG /* returns x[i] */ #define v_entry(x,i) (((i) < 0 || (i) >= (x)->dim) ? \ error(E_BOUNDS,"v_entry"), 0.0 : (x)->ve[i] ) /* x[i] <- val */ #define v_set_val(x,i,val) ((x)->ve[i] = ((i) < 0 || (i) >= (x)->dim) ? \ error(E_BOUNDS,"v_set_val"), 0.0 : (val)) /* x[i] <- x[i] + val */ #define v_add_val(x,i,val) ((x)->ve[i] += ((i) < 0 || (i) >= (x)->dim) ? \ error(E_BOUNDS,"v_add_val"), 0.0 : (val)) /* x[i] <- x[i] - val */ #define v_sub_val(x,i,val) ((x)->ve[i] -= ((i) < 0 || (i) >= (x)->dim) ? \ error(E_BOUNDS,"v_sub_val"), 0.0 : (val)) /* returns A[i][j] */ #define m_entry(A,i,j) (((i) < 0 || (i) >= (A)->m || \ (j) < 0 || (j) >= (A)->n) ? \ error(E_BOUNDS,"m_entry"), 0.0 : (A)->me[i][j] ) /* A[i][j] <- val */ #define m_set_val(A,i,j,val) ((A)->me[i][j] = ((i) < 0 || (i) >= (A)->m || \ (j) < 0 || (j) >= (A)->n) ? \ error(E_BOUNDS,"m_set_val"), 0.0 : (val) ) /* A[i][j] <- A[i][j] + val */ #define m_add_val(A,i,j,val) ((A)->me[i][j] += ((i) < 0 || (i) >= (A)->m || \ (j) < 0 || (j) >= (A)->n) ? \ error(E_BOUNDS,"m_add_val"), 0.0 : (val) ) /* A[i][j] <- A[i][j] - val */ #define m_sub_val(A,i,j,val) ((A)->me[i][j] -= ((i) < 0 || (i) >= (A)->m || \ (j) < 0 || (j) >= (A)->n) ? \ error(E_BOUNDS,"m_sub_val"), 0.0 : (val) ) #else /* returns x[i] */ #define v_entry(x,i) ((x)->ve[i]) /* x[i] <- val */ #define v_set_val(x,i,val) ((x)->ve[i] = (val)) /* x[i] <- x[i] + val */ #define v_add_val(x,i,val) ((x)->ve[i] += (val)) /* x[i] <- x[i] - val */ #define v_sub_val(x,i,val) ((x)->ve[i] -= (val)) /* returns A[i][j] */ #define m_entry(A,i,j) ((A)->me[i][j]) /* A[i][j] <- val */ #define m_set_val(A,i,j,val) ((A)->me[i][j] = (val) ) /* A[i][j] <- A[i][j] + val */ #define m_add_val(A,i,j,val) ((A)->me[i][j] += (val) ) /* A[i][j] <- A[i][j] - val */ #define m_sub_val(A,i,j,val) ((A)->me[i][j] -= (val) ) #endif /* I/O routines */ #ifndef ANSI_C extern void v_foutput(),m_foutput(),px_foutput(); extern void iv_foutput(); extern VEC *v_finput(); extern MAT *m_finput(); extern PERM *px_finput(); extern IVEC *iv_finput(); extern int fy_or_n(), fin_int(), yn_dflt(), skipjunk(); extern double fin_double(); #else /* print x on file fp */ void v_foutput(FILE *fp,VEC *x), /* print A on file fp */ m_foutput(FILE *fp,MAT *A), /* print px on file fp */ px_foutput(FILE *fp,PERM *px); /* print ix on file fp */ void iv_foutput(FILE *fp,IVEC *ix); /* Note: if out is NULL, then returned object is newly allocated; Also: if out is not NULL, then that size is assumed */ /* read in vector from fp */ VEC *v_finput(FILE *fp,VEC *out); /* read in matrix from fp */ MAT *m_finput(FILE *fp,MAT *out); /* read in permutation from fp */ PERM *px_finput(FILE *fp,PERM *out); /* read in int vector from fp */ IVEC *iv_finput(FILE *fp,IVEC *out); /* fy_or_n -- yes-or-no to question in string s -- question written to stderr, input from fp -- if fp is NOT a tty then return y_n_dflt */ int fy_or_n(FILE *fp,char *s); /* yn_dflt -- sets the value of y_n_dflt to val */ int yn_dflt(int val); /* fin_int -- return integer read from file/stream fp -- prompt s on stderr if fp is a tty -- check that x lies between low and high: re-prompt if fp is a tty, error exit otherwise -- ignore check if low > high */ int fin_int(FILE *fp,char *s,int low,int high); /* fin_double -- return double read from file/stream fp -- prompt s on stderr if fp is a tty -- check that x lies between low and high: re-prompt if fp is a tty, error exit otherwise -- ignore check if low > high */ double fin_double(FILE *fp,char *s,double low,double high); /* it skips white spaces and strings of the form #....\n Here .... is a comment string */ int skipjunk(FILE *fp); #endif /* MACROS */ /* macros to use stdout and stdin instead of explicit fp */ #define v_output(vec) v_foutput(stdout,vec) #define v_input(vec) v_finput(stdin,vec) #define m_output(mat) m_foutput(stdout,mat) #define m_input(mat) m_finput(stdin,mat) #define px_output(px) px_foutput(stdout,px) #define px_input(px) px_finput(stdin,px) #define iv_output(iv) iv_foutput(stdout,iv) #define iv_input(iv) iv_finput(stdin,iv) /* general purpose input routine; skips comments # ... \n */ #define finput(fp,prompt,fmt,var) \ ( ( isatty(fileno(fp)) ? fprintf(stderr,prompt) : skipjunk(fp) ), \ fscanf(fp,fmt,var) ) #define input(prompt,fmt,var) finput(stdin,prompt,fmt,var) #define fprompter(fp,prompt) \ ( isatty(fileno(fp)) ? fprintf(stderr,prompt) : skipjunk(fp) ) #define prompter(prompt) fprompter(stdin,prompt) #define y_or_n(s) fy_or_n(stdin,s) #define in_int(s,lo,hi) fin_int(stdin,s,lo,hi) #define in_double(s,lo,hi) fin_double(stdin,s,lo,hi) /* Copying routines */ #ifndef ANSI_C extern MAT *_m_copy(), *m_move(), *vm_move(); extern VEC *_v_copy(), *v_move(), *mv_move(); extern PERM *px_copy(); extern IVEC *iv_copy(), *iv_move(); extern BAND *bd_copy(); #else /* copy in to out starting at out[i0][j0] */ extern MAT *_m_copy(MAT *in,MAT *out,int i0,int j0), * m_move(MAT *in, int, int, int, int, MAT *out, int, int), *vm_move(VEC *in, int, MAT *out, int, int, int, int); /* copy in to out starting at out[i0] */ extern VEC *_v_copy(VEC *in,VEC *out,int i0), * v_move(VEC *in, int, int, VEC *out, int), *mv_move(MAT *in, int, int, int, int, VEC *out, int); extern PERM *px_copy(PERM *in,PERM *out); extern IVEC *iv_copy(IVEC *in,IVEC *out), *iv_move(IVEC *in, int, int, IVEC *out, int); extern BAND *bd_copy(BAND *in,BAND *out); #endif /* MACROS */ #define m_copy(in,out) _m_copy(in,out,0,0) #define v_copy(in,out) _v_copy(in,out,0) /* function prototypes */ MAT *im_finput( FILE* fp, MAT* a ); MAT *bm_finput( FILE* fp, MAT* a ); PERM *ipx_finput(FILE* fp, PERM* px); PERM *bpx_finput(FILE* fp, PERM* px); VEC *ifin_vec( FILE* fp, VEC* x ); VEC *bfin_vec( FILE* fp, VEC* x ); /* Initialisation routines -- to be zero, ones, random or identity */ VEC *v_zero(VEC *); VEC *v_rand(VEC *); VEC *v_ones(VEC *); MAT *m_zero(MAT *), *m_ident(MAT *), *m_rand(MAT *), *m_ones(MAT *); PERM *px_ident(PERM *); IVEC *iv_zero(IVEC *); /* Basic vector operations */ #ifndef ANSI_C extern VEC *sv_mlt(), *mv_mlt(), *vm_mlt(), *v_add(), *v_sub(), *px_vec(), *pxinv_vec(), *v_mltadd(), *v_map(), *_v_map(), *v_lincomb(), *v_linlist(); extern double v_min(), v_max(), v_sum(); extern VEC *v_star(), *v_slash(), *v_sort(); extern double _in_prod(), __ip__(); extern void __mltadd__(), __add__(), __sub__(), __smlt__(), __zero__(); #else extern VEC *sv_mlt(double,VEC *,VEC *), /* out <- s.x */ *mv_mlt(MAT *,VEC *,VEC *), /* out <- A.x */ *vm_mlt(MAT *,VEC *,VEC *), /* out^T <- x^T.A */ *v_add(VEC *,VEC *,VEC *), /* out <- x + y */ *v_sub(VEC *,VEC *,VEC *), /* out <- x - y */ *px_vec(PERM *,VEC *,VEC *), /* out <- P.x */ *pxinv_vec(PERM *,VEC *,VEC *), /* out <- P^{-1}.x */ *v_mltadd(VEC *,VEC *,double,VEC *), /* out <- x + s.y */ #ifdef PROTOTYPES_IN_STRUCT *v_map(double (*f)(double),VEC *,VEC *), /* out[i] <- f(x[i]) */ *_v_map(double (*f)(void *,double),void *,VEC *,VEC *), #else *v_map(double (*f)(),VEC *,VEC *), /* out[i] <- f(x[i]) */ *_v_map(double (*f)(),void *,VEC *,VEC *), #endif *v_lincomb(int,VEC **,Real *,VEC *), /* out <- sum_i s[i].x[i] */ *v_linlist(VEC *out,VEC *v1,double a1,...); /* out <- s1.x1 + s2.x2 + ... */ /* returns min_j x[j] (== x[i]) */ extern double v_min(VEC *, int *), /* returns max_j x[j] (== x[i]) */ v_max(VEC *, int *), /* returns sum_i x[i] */ v_sum(VEC *); /* Hadamard product: out[i] <- x[i].y[i] */ extern VEC *v_star(VEC *, VEC *, VEC *), /* out[i] <- x[i] / y[i] */ *v_slash(VEC *, VEC *, VEC *), /* sorts x, and sets order so that sorted x[i] = x[order[i]] */ *v_sort(VEC *, PERM *); /* returns inner product starting at component i0 */ extern double _in_prod(VEC *x,VEC *y,int i0), /* returns sum_{i=0}^{len-1} x[i].y[i] */ __ip__(Real *,Real *,int); /* see v_mltadd(), v_add(), v_sub() and v_zero() */ void __mltadd__(Real *,Real *,double,int), __add__(Real *,Real *,Real *,int), __sub__(Real *,Real *,Real *,int), __smlt__(Real *,double,Real *,int), __zero__(Real *,int); #endif /* MACRO */ /* usual way of computing the inner product */ #define in_prod(a,b) _in_prod(a,b,0) /* Norms */ /* scaled vector norms -- scale == NULL implies unscaled */ #ifndef ANSI_C extern double _v_norm1(), _v_norm2(), _v_norm_inf(), m_norm1(), m_norm_inf(), m_norm_frob(); #else /* returns sum_i |x[i]/scale[i]| */ extern double _v_norm1(VEC *x,VEC *scale), /* returns (scaled) Euclidean norm */ _v_norm2(VEC *x,VEC *scale), /* returns max_i |x[i]/scale[i]| */ _v_norm_inf(VEC *x,VEC *scale); /* unscaled matrix norms */ extern double m_norm1(MAT *A), m_norm_inf(MAT *A), m_norm_frob(MAT *A); #endif /* MACROS */ /* unscaled vector norms */ #define v_norm1(x) _v_norm1(x,VNULL) #define v_norm2(x) _v_norm2(x,VNULL) #define v_norm_inf(x) _v_norm_inf(x,VNULL) /* Basic matrix operations */ #ifndef ANSI_C extern MAT *sm_mlt(), *m_mlt(), *mmtr_mlt(), *mtrm_mlt(), *m_add(), *m_sub(), *sub_mat(), *m_transp(), *ms_mltadd(); extern BAND *bd_transp(); extern MAT *px_rows(), *px_cols(), *swap_rows(), *swap_cols(), *_set_row(), *_set_col(); extern VEC *get_row(), *get_col(), *sub_vec(), *mv_mltadd(), *vm_mltadd(); #else extern MAT *sm_mlt(double s,MAT *A,MAT *out), /* out <- s.A */ *m_mlt(MAT *A,MAT *B,MAT *out), /* out <- A.B */ *mmtr_mlt(MAT *A,MAT *B,MAT *out), /* out <- A.B^T */ *mtrm_mlt(MAT *A,MAT *B,MAT *out), /* out <- A^T.B */ *m_add(MAT *A,MAT *B,MAT *out), /* out <- A + B */ *m_sub(MAT *A,MAT *B,MAT *out), /* out <- A - B */ *sub_mat(MAT *A,int,int,int,int,MAT *out), *m_transp(MAT *A,MAT *out), /* out <- A^T */ /* out <- A + s.B */ *ms_mltadd(MAT *A,MAT *B,double s,MAT *out); extern BAND *bd_transp(BAND *in, BAND *out); /* out <- A^T */ extern MAT *px_rows(PERM *px,MAT *A,MAT *out), /* out <- P.A */ *px_cols(PERM *px,MAT *A,MAT *out), /* out <- A.P^T */ *swap_rows(MAT *,int,int,int,int), *swap_cols(MAT *,int,int,int,int), /* A[i][j] <- out[j], j >= j0 */ *_set_col(MAT *A,int i,VEC *out,int j0), /* A[i][j] <- out[i], i >= i0 */ *_set_row(MAT *A,int j,VEC *out,int i0); extern VEC *get_row(MAT *,int,VEC *), *get_col(MAT *,int,VEC *), *sub_vec(VEC *,int,int,VEC *), /* out <- x + s.A.y */ *mv_mltadd(VEC *x,VEC *y,MAT *A,double s,VEC *out), /* out^T <- x^T + s.y^T.A */ *vm_mltadd(VEC *x,VEC *y,MAT *A,double s,VEC *out); #endif /* MACROS */ /* row i of A <- vec */ #define set_row(mat,row,vec) _set_row(mat,row,vec,0) /* col j of A <- vec */ #define set_col(mat,col,vec) _set_col(mat,col,vec,0) /* Basic permutation operations */ #ifndef ANSI_C extern PERM *px_mlt(), *px_inv(), *px_transp(); extern int px_sign(); #else extern PERM *px_mlt(PERM *px1,PERM *px2,PERM *out), /* out <- px1.px2 */ *px_inv(PERM *px,PERM *out), /* out <- px^{-1} */ /* swap px[i] and px[j] */ *px_transp(PERM *px,int i,int j); /* returns sign(px) = +1 if px product of even # transpositions -1 if ps product of odd # transpositions */ extern int px_sign(PERM *); #endif /* Basic integer vector operations */ #ifndef ANSI_C extern IVEC *iv_add(), *iv_sub(), *iv_sort(); #else extern IVEC *iv_add(IVEC *ix,IVEC *iy,IVEC *out), /* out <- ix + iy */ *iv_sub(IVEC *ix,IVEC *iy,IVEC *out), /* out <- ix - iy */ /* sorts ix & sets order so that sorted ix[i] = old ix[order[i]] */ *iv_sort(IVEC *ix, PERM *order); #endif /* miscellaneous functions */ #ifndef ANSI_C extern double square(), cube(), mrand(); extern void smrand(), mrandlist(); extern void m_dump(), px_dump(), v_dump(), iv_dump(); extern MAT *band2mat(); extern BAND *mat2band(); #else double square(double x), /* returns x^2 */ cube(double x), /* returns x^3 */ mrand(void); /* returns random # in [0,1) */ void smrand(int seed), /* seeds mrand() */ mrandlist(Real *x, int len); /* generates len random numbers */ void m_dump(FILE *fp,MAT *a), px_dump(FILE *,PERM *px), v_dump(FILE *fp,VEC *x), iv_dump(FILE *fp, IVEC *ix); MAT *band2mat(BAND *bA, MAT *A); BAND *mat2band(MAT *A, int lb,int ub, BAND *bA); #endif /* miscellaneous constants */ #define VNULL ((VEC *)NULL) #define MNULL ((MAT *)NULL) #define PNULL ((PERM *)NULL) #define IVNULL ((IVEC *)NULL) #define BDNULL ((BAND *)NULL) /* varying number of arguments */ #ifdef ANSI_C #include /* prototypes */ int v_get_vars(int dim,...); int iv_get_vars(int dim,...); int m_get_vars(int m,int n,...); int px_get_vars(int dim,...); int v_resize_vars(int new_dim,...); int iv_resize_vars(int new_dim,...); int m_resize_vars(int m,int n,...); int px_resize_vars(int new_dim,...); int v_free_vars(VEC **,...); int iv_free_vars(IVEC **,...); int px_free_vars(PERM **,...); int m_free_vars(MAT **,...); #elif VARARGS /* old varargs is used */ #include /* prototypes */ int v_get_vars(); int iv_get_vars(); int m_get_vars(); int px_get_vars(); int v_resize_vars(); int iv_resize_vars(); int m_resize_vars(); int px_resize_vars(); int v_free_vars(); int iv_free_vars(); int px_free_vars(); int m_free_vars(); #endif #endif /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* Header file for ``matrix2.a'' library file */ #ifndef MATRIX2H #define MATRIX2H /* Unless otherwise specified, factorisation routines overwrite the matrix that is being factorised */ #ifndef ANSI_C extern MAT *BKPfactor(), *CHfactor(), *LUfactor(), *QRfactor(), *QRCPfactor(), *LDLfactor(), *Hfactor(), *MCHfactor(), *m_inverse(); extern double LUcondest(), QRcondest(); extern MAT *makeQ(), *makeR(), *makeHQ(), *makeH(); extern MAT *LDLupdate(), *QRupdate(); extern VEC *BKPsolve(), *CHsolve(), *LUsolve(), *_Qsolve(), *QRsolve(), *LDLsolve(), *Usolve(), *Lsolve(), *Dsolve(), *LTsolve(), *UTsolve(), *LUTsolve(), *QRCPsolve(); extern BAND *bdLUfactor(), *bdLDLfactor(); extern VEC *bdLUsolve(), *bdLDLsolve(); extern VEC *hhvec(); extern VEC *hhtrvec(); extern MAT *hhtrrows(); extern MAT *hhtrcols(); extern void givens(); extern VEC *rot_vec(); /* in situ */ extern MAT *rot_rows(); /* in situ */ extern MAT *rot_cols(); /* in situ */ /* eigenvalue routines */ extern VEC *trieig(), *symmeig(); extern MAT *schur(); extern void schur_evals(); extern MAT *schur_vecs(); /* singular value decomposition */ extern VEC *bisvd(), *svd(); /* matrix powers and exponent */ MAT *_m_pow(); MAT *m_pow(); MAT *m_exp(), *_m_exp(); MAT *m_poly(); /* FFT */ void fft(); void ifft(); #else /* forms Bunch-Kaufman-Parlett factorisation for symmetric indefinite matrices */ extern MAT *BKPfactor(MAT *A,PERM *pivot,PERM *blocks), /* Cholesky factorisation of A (symmetric, positive definite) */ *CHfactor(MAT *A), /* LU factorisation of A (with partial pivoting) */ *LUfactor(MAT *A,PERM *pivot), /* QR factorisation of A; need dim(diag) >= # rows of A */ *QRfactor(MAT *A,VEC *diag), /* QR factorisation of A with column pivoting */ *QRCPfactor(MAT *A,VEC *diag,PERM *pivot), /* L.D.L^T factorisation of A */ *LDLfactor(MAT *A), /* Hessenberg factorisation of A -- for schur() */ *Hfactor(MAT *A,VEC *diag1,VEC *diag2), /* modified Cholesky factorisation of A; actually factors A+D, D diagonal with no diagonal entry in the factor < sqrt(tol) */ *MCHfactor(MAT *A,double tol), *m_inverse(MAT *A,MAT *out); /* returns condition estimate for A after LUfactor() */ extern double LUcondest(MAT *A,PERM *pivot), /* returns condition estimate for Q after QRfactor() */ QRcondest(MAT *A); /* Note: The make..() and ..update() routines assume that the factorisation has already been carried out */ /* Qout is the "Q" (orthongonal) matrix from QR factorisation */ extern MAT *makeQ(MAT *A,VEC *diag,MAT *Qout), /* Rout is the "R" (upper triangular) matrix from QR factorisation */ *makeR(MAT *A,MAT *Rout), /* Qout is orthogonal matrix in Hessenberg factorisation */ *makeHQ(MAT *A,VEC *diag1,VEC *diag2,MAT *Qout), /* Hout is the Hessenberg matrix in Hessenberg factorisation */ *makeH(MAT *A,MAT *Hout); /* updates L.D.L^T factorisation for A <- A + alpha.u.u^T */ extern MAT *LDLupdate(MAT *A,VEC *u,double alpha), /* updates QR factorisation for QR <- Q.(R+u.v^T) Note: we need explicit Q & R matrices, from makeQ() and makeR() */ *QRupdate(MAT *Q,MAT *R,VEC *u,VEC *v); /* Solve routines assume that the corresponding factorisation routine has already been applied to the matrix along with auxiliary objects (such as pivot permutations) These solve the system A.x = b, except for LUTsolve and QRTsolve which solve the transposed system A^T.x. = b. If x is NULL on entry, then it is created. */ extern VEC *BKPsolve(MAT *A,PERM *pivot,PERM *blocks,VEC *b,VEC *x), *CHsolve(MAT *A,VEC *b,VEC *x), *LDLsolve(MAT *A,VEC *b,VEC *x), *LUsolve(MAT *A,PERM *pivot,VEC *b,VEC *x), *_Qsolve(MAT *A,VEC *,VEC *,VEC *, VEC *), *QRsolve(MAT *A,VEC *,VEC *b,VEC *x), *QRTsolve(MAT *A,VEC *,VEC *b,VEC *x), /* Triangular equations solve routines; U for upper triangular, L for lower traingular, D for diagonal if diag_val == 0.0 use that values in the matrix */ *Usolve(MAT *A,VEC *b,VEC *x,double diag_val), *Lsolve(MAT *A,VEC *b,VEC *x,double diag_val), *Dsolve(MAT *A,VEC *b,VEC *x), *LTsolve(MAT *A,VEC *b,VEC *x,double diag_val), *UTsolve(MAT *A,VEC *b,VEC *x,double diag_val), *LUTsolve(MAT *A,PERM *,VEC *,VEC *), *QRCPsolve(MAT *QR,VEC *diag,PERM *pivot,VEC *b,VEC *x); extern BAND *bdLUfactor(BAND *A,PERM *pivot), *bdLDLfactor(BAND *A); extern VEC *bdLUsolve(BAND *A,PERM *pivot,VEC *b,VEC *x), *bdLDLsolve(BAND *A,VEC *b,VEC *x); extern VEC *hhvec(VEC *,int,Real *,VEC *,Real *); extern VEC *hhtrvec(VEC *,double,int,VEC *,VEC *); extern MAT *hhtrrows(MAT *,int,int,VEC *,double); extern MAT *hhtrcols(MAT *,int,int,VEC *,double); extern void givens(double,double,Real *,Real *); extern VEC *rot_vec(VEC *,int,int,double,double,VEC *); /* in situ */ extern MAT *rot_rows(MAT *,int,int,double,double,MAT *); /* in situ */ extern MAT *rot_cols(MAT *,int,int,double,double,MAT *); /* in situ */ /* eigenvalue routines */ /* compute eigenvalues of tridiagonal matrix with diagonal entries a[i], super & sub diagonal entries b[i]; eigenvectors stored in Q (if not NULL) */ extern VEC *trieig(VEC *a,VEC *b,MAT *Q), /* sets out to be vector of eigenvectors; eigenvectors stored in Q (if not NULL). A is unchanged */ *symmeig(MAT *A,MAT *Q,VEC *out); /* computes real Schur form = Q^T.A.Q */ extern MAT *schur(MAT *A,MAT *Q); /* computes real and imaginary parts of the eigenvalues of A after schur() */ extern void schur_evals(MAT *A,VEC *re_part,VEC *im_part); /* computes real and imaginary parts of the eigenvectors of A after schur() */ extern MAT *schur_vecs(MAT *T,MAT *Q,MAT *X_re,MAT *X_im); /* singular value decomposition */ /* computes singular values of bi-diagonal matrix with diagonal entries a[i] and superdiagonal entries b[i]; singular vectors stored in U and V (if not NULL) */ VEC *bisvd(VEC *a,VEC *b,MAT *U,MAT *V), /* sets out to be vector of singular values; singular vectors stored in U and V */ *svd(MAT *A,MAT *U,MAT *V,VEC *out); /* matrix powers and exponent */ MAT *_m_pow(MAT *,int,MAT *,MAT *); MAT *m_pow(MAT *,int, MAT *); MAT *m_exp(MAT *,double,MAT *); MAT *_m_exp(MAT *,double,MAT *,int *,int *); MAT *m_poly(MAT *,VEC *,MAT *); /* FFT */ void fft(VEC *,VEC *); void ifft(VEC *,VEC *); #endif #endif /* From err.h */ /* error(E_TYPE,"myfunc") raises error type E_TYPE for function my_func() */ #define error(err_num,fn_name) ev_err(__FILE__,err_num,__LINE__,fn_name,0) /* End of Header files *****************************************************/ int Optopt; extern int errno; void usage(void){printf(pstUsage);} void die(char* pst){ printf("%s, with error#%d in testcfg\n",pst,errno); exit(1); } int ev_err(char *stFile,int Num,int Line,char * Func,int Type); int readOptions(int argc, char*argv[]); /* interpret command line */ int checkWriteCfg(char *pstIni); /* test and/or write the config file */ int doEdit(char *pstIni); /* edit and/or write the config file */ int readProcConfig(char* pstIni); /* read and process the config file */ double ran1(long *pId) ;/* uniform variates */ double sran1(long *pId); /* normal variates */ void copyarrvec(float* from, VEC * result, int max); void copyarrmat(float* from, MAT * result, int maxR, int maxC); int symsqrtmatrix(MAT *mS, MAT *mC, MAT * mCheck); void addvec(VEC * result, VEC * add); /* adds add to result */ void printarr(float * vr, int Max); /* prints a vector */ void printvec(VEC * vector); void printmat(MAT * matrix); enum { vMEAN,vSHIFT,mSIGMA,mPHI,N_OBS,N_IC,SEED }; struct CfgStrings Cfg[] = { { "mean" , NULL, NULL }, { "shift" , NULL, NULL }, { "sigma" , NULL, NULL }, { "phi" , NULL, NULL }, { "n_obs" , NULL, NULL }, { "n_incontrol" , NULL, NULL }, { "seed" , NULL, NULL }, { NULL , NULL, NULL } /* array terminator. REQUIRED !!! */ }; /* options and globals: */ const char * pstRCSID="# $Id: var1.c,v 1.20 2000/02/01 21:11:27 drf5n Exp drf5n $\n\ "; short fDoEdit = 0; short fVerbose = 0; short fTerse = 0; char *pst, *pstIni; long int maxIC=10; /* upper bound on in control observations */ long int maxObs=20; /* upper bound on number of total observations */ int scDim; VEC *vcMean = NULL; /* process mean */ VEC *vcShift = NULL; /* process mean */ MAT *mSigma = NULL; /* Noise Covariance Meschach form */ MAT *mC = NULL; /* square root Noise Covariance Meschach form */ MAT *mPhi = NULL; /* process VAR autocorrelation parameter */ long int seed = 0; /* random number seed */ MAT *mXXX = NULL; /* testing */ /* the program *********************************************************/ int main(int argc, char*argv[]){ short fInCtrl = 1; int scRow,scCol; int i; size_t scLen; long scObs; VEC *vcN; /* Random vector N(0,mSigma) */ VEC *vcE; /* Random vector N(0,I) */ VEC *vcX; /* Random vector N(vcMean,mSigma*(I-mPhi) */ VEC *vcXpast; /* Random vector N(vcMean,mSigma*(I-mPhi) */ VEC *vcY; /* Random vector N(vcMean,mSigma*(I-mPhi) */ scRow = 1; scCol = 1; scDim = 0; scObs = 0; pstIni="var1.ini"; pst = NULL; scLen =0; if (readOptions(argc, argv) != 0 ) /* interpret command line */ die("Options not understood"); if(checkWriteCfg(pstIni) != 0) /* test and/or write the config file */ die("Couldn't read or write the configuration file"); if (fDoEdit) if (doEdit( pstIni) !=0 ) /* edit and/or write the config file */ printf("error during edit of the configuration file\n"); if(readProcConfig( pstIni) !=0) /* read and process the config file */ die("Couldn't read configuration file"); vcN = v_get(scDim); vcE = v_get(scDim); vcX = v_get(scDim); vcXpast = v_get(scDim); vcY = v_get(scDim); /* the configuration is now fixed */ if(fVerbose){ printf("#begin VERBOSE \n"); printf("%s \n",pstRCSID); printf("initialization file = %s",pstIni); printf("\nobservations = %ld",maxObs); printf("\nnumber in control = %ld ",maxIC); printf("\nseed = %ld ",-seed); printf("\nmean = ");printvec(vcMean); printf("\n\nshift = ");printvec(vcShift); printf("\n\ncovariance = \"");printmat(mSigma); printf("\"\n\nautocorrelation = \"");printmat(mPhi); } /* setup for calcs */ if (symsqrtmatrix(mSigma, mC, mXXX) != 0) { for(i=0;ime[1][i]); printf("= eigenvalues\n\n"); printmat(mC);printf(" ^== eigenvectors\n"); die("Error: The symmetric square root of the noise covariance\n\ does not exist\n"); } #if DEBUG > 9 printmat(mXXX); printf(" --- Check -CC+S=0 ??? \n"); #endif if(fVerbose){ printf("\n\nC = \"");printmat(mC); printf("\"\n\n"); printf("#end VERBOSE\n"); } v_zero(vcX); /* initialize the VAR() series */ for (scObs=0;scObsve[i]=sran1(&seed); /* vcN = N(0,I) */ } if(!fInCtrl) addvec(vcN,vcShift); /* vcN = N(vcShift,I) */ mv_mlt(mC,vcN,vcE); /* e = CN */ mv_mlt(mPhi,vcXpast,vcX); /* X_t = Phi X_{t-1} */ addvec(vcX,vcE); /* X_t = Phi X_{t-1} + e */ vcY=v_copy(vcX,vcY); addvec(vcY,vcMean); /* Y_t = X_t +u */ printvec(vcY); /* process */ printf("\n"); } return 0; } /* end of main ******************************/ /* *************** drf routines ************* */ void printvec(VEC * v){ int i; if (v->dim > 1 ) printf("%f",v->ve[0]); for(i=1;i< v->dim;++i) printf(" %f",v->ve[i]); } void printmat(MAT * m){ int iR,iC; for(iR=0;iR< m->m;iR++) { if(iR != 0 ) printf("\n"); for(iC=0;iCn;iC++) printf("%f ",m->me[iR][iC]); } } void printarr(float arr[], int max){ int i; for(i=0;idim;i++) result->ve[i]+=add->ve[i]; } void copyarrvec(float* from, VEC * result, int max){ /* copies a x[0...n-1] c array to a Meschach Vector */ int i; #if DEBUG > 9 printf("in copyarrvec"); printf("from[0]: %f [1]: %f",from[0],from[1]); #endif for(i=0;ive[i]=from[i]; } void copyarrmat(float* from, MAT * result, int maxR, int maxC){ /* copies a x[0...n-1] c array to a Meschach RxC array * blindly copies, whether there is data there or not */ int iR, iC; #if DEBUG > 9 printf("in copyarrmat"); #endif for(iR=0;iRme[iR][iC]=from[iR*maxC+iC]; } void matscalarMultiply(MAT * result, double scalar, int maxR, int maxC){ int iR, iC; #if DEBUG > 9 printf("in matscalarMult"); #endif for(iR=0;iRme[iR][iC] *= scalar; } void matmatAdd(MAT * result, MAT * from, int maxR, int maxC){ int iR, iC; #if DEBUG > 9 printf("in matmatAdd"); #endif for(iR=0;iRme[iR][iC]+=from->me[iR][iC]; } void matmatMultiply(MAT * result, MAT *A, MAT * B, int maxR, int maxI,int maxC) { int iR,iI, iC; #if DEBUG > 9 printf("in matmatMultiply"); #endif for(iR=0;iRme[iR][iC]=0.0; for(iI=0;iIme[iR][iC] += A->me[iR][iI] * B->me[iI][iC]; ; } } int symsqrtmatrix(MAT *mS, MAT *mC, MAT * mCheck){ /* Find C such that CC=S by eigenvector/value decomp * I = a check: S-CC * by: * S = V'AV * where A = eigenvals V = eigenvecs * where the square roots may be taken on each eigenval * C= V(D^0.5)V' * so CC = S * * First, the matlab code: * %find the symmetric square root of mSig (Rencher 1996 p 412) * [V D] = eig(mSig); mSigH = (V * ( D ^.5))*V' ; */ int i, iR, iC, n; MAT * mW; MAT * mX; MAT * mV; VEC * vA; VEC * vD; n = mS->n; if (mS->n != mS->m) { perror ("symsqrtmatrix: matrix not square"); return -1; } mW=m_get(n,n); /* temporary */ mX=m_get(n,n); /* temporary */ mV=m_get(n,n); /* eigenvectors */ vA=v_get(n); /* eigenvalues */ vD=v_get(n); /* eigenvalues^.5 */ m_copy(mS,mW); symmeig(mS,mV,vA); #if DEBUG > 9 printf("in symsqrtmatrix():\n"); printvec(vA );printf(" <-- eigenvalues\n"); printmat(mV);printf(" ^-- eigenvectors\n"); #endif m_zero(mW); for(i=0;ive[i]<0) { /* failure */ mC=m_copy(mV,mC); for(iC=0;iCme[1][iC]=vA->ve[iC]; return -1; } vD->ve[i]=sqrt(vA->ve[i]); /* square root eigenvalues */ mW->me[i][i]=vD->ve[i]; } matmatMultiply(mX,mV,mW,n,n,n); /* = V*D(^.5) */ for(iR=0; iRme[iR][iC]=mV->me[iC][iR] ; /* transpose */ } } matmatMultiply(mC,mX,mW,n,n,n); /* (V * ( D ^.5))*V' */ /* check */ matmatMultiply(mCheck,mC,mC,n,n,n); #if DEBUG > 9 printmat(mCheck);printf(":CC\n"); /* drf5n */ #endif matscalarMultiply(mCheck,-1.0,n,n); matmatAdd(mCheck,mS,n,n); /* should be zero */ #if DEBUG > 9 printmat(mS);printf(" ^-- covariance in lower\n"); printvec(vA);printf(" <-- eigenvalues\n"); printmat(mV);printf(" ^-- eigenvectors\n"); printmat(mC); printf(" ^-- sqrt(covar)\n"); printmat(mCheck); printf(" ^-- Check Matrix - should = 0\n"); #endif return 0; } /**************** RANDOM NUMBERS **************/ /* uniform random variate */ double ran1(long *pId) {/* */ long k1; k1 = *pId/127773; if(*pId == 0 ) *pId=1; *pId=16807 * ( *pId -k1*127773) -k1 * 2836; if( *pId < 0 ) *pId += 2147483647; #if DEBUG > 9 printf("ran1: %f\n",((double)*pId)*(4.656612875e-10)); #endif return ((double)*pId)*(4.656612875e-10); } /* standard normal random variate */ double sran1(long *pId) {/* */ static int fExtra=0; static double X2; double fac,rsq,v1,v2; if (pId == NULL) {perror("Null pointer in sran1");exit(1);} if (*pId < 0) fExtra=0; if (fExtra == 0) { do { v1=2.0*ran1(pId)-1.0; v2=2.0*ran1(pId)-1.0; rsq=v1*v1+v2*v2; #if DEBUG > 9 printf("looking for standard normals %f\n",rsq); #endif } while (rsq >= 1.0 || rsq == 0.0); fac=sqrt(-2.0*log(rsq)/rsq); X2=v1*fac; fExtra=1; return v2*fac; } else { fExtra=0; return X2; } } int readOptions(int argc, char*argv[]){ int c; while ( (c=getopt(argc,argv,"ti:eh?n:k:r:s:v")) != EOF && c != 0) switch (c) { case 'i': pstIni = optarg; break; case 'v': fVerbose = 1; break; case 't': fTerse = 1; break; case 'e': fDoEdit = 1; break; case 'k': /* number of in control observations */ maxIC = atol(optarg); Cfg[N_IC].def = optarg; break; case 'r': case 'n': /* number of observations */ maxObs = atol(optarg); Cfg[N_OBS].def = optarg; if (maxObs <1) { printf("Need more than n=%li observations\n", maxObs); return 1; } break; case 's': /* random number seed */ seed = atol(optarg); Cfg[SEED].def = optarg; break; case '?': case 'h': usage(); return 1; default: if (isprint (Optopt)) fprintf (stderr, "Unknown option `-%c'.\n", Optopt); else fprintf (stderr, "Unknown option character `\\x%x'.\n", Optopt); return 1; } return 0; } int checkWriteCfg(char *pstIni){ /* test and/or write the config file */ FILE * pf; if((pf =fopen(pstIni,"r")) !=NULL) fclose (pf); else { if((pf =fopen(pstIni,"w")) == NULL) { printf("couldn't open %s\n",pstIni); return 1; } else { fputs(pstIniDefault,pf); if (fclose(pf) != 0){ die("error closing config file"); } printf("# Wrote default configuration file to %s\n",pstIni); } } return 0; } int doEdit(char *pstIni){ /* edit and/or write the config file */ char stBuffer[1024]; int retval; if(getenv("EDITOR")!=NULL){ sprintf(stBuffer,"%s %s",getenv("EDITOR"), pstIni); }else{ printf("The environmental variable EDITOR not set. Please set it to your\n\ preferred text editor. \nAttempting to use \"%s\"\n",DEFAULT_EDITOR); sprintf(stBuffer,"%s %s",DEFAULT_EDITOR, pstIni); } retval=(system(stBuffer)); if(retval != 0 ){ printf("Editing: \"%s\" gave exit code %d\n",stBuffer,retval); } return retval; } int readProcConfig(char* pstIni){ /* read and process the config file */ int sc, maxarr = 0; long cX; float *arr = NULL; /* array */ if((cX=CfgRead(pstIni , Cfg)) == 0) return 1; else { /* The command line args take precedence over the config file */ if(Cfg[N_OBS].def) { maxObs = atol(Cfg[N_OBS].def );} else if(Cfg[N_OBS].data) { maxObs = atol(Cfg[N_OBS].data);} if(Cfg[N_IC].def) { maxIC = atol(Cfg[N_IC].def );} else if(Cfg[N_IC].data) { maxIC = atol(Cfg[N_IC].data);} if(Cfg[SEED].def) { seed = atol(Cfg[SEED].def );} else if(Cfg[SEED].data) { seed = atol(Cfg[SEED].data);} seed = -labs(seed); /* make the seed meaningful to ran1() */ scDim = readfloats( Cfg[vMEAN].data, &arr,&maxarr); if (scDim <= 0 ) { printf("Mean Dimension (%s) too small\n",Cfg[vMEAN].data); die(""); } #if DEBUG > 9 printf("read mean: "); printarr(arr,scDim); #endif /* get some memory */ vcMean = v_get(scDim); vcShift = v_get(scDim); mSigma = m_get(scDim,scDim); mC = m_get(scDim,scDim); mPhi = m_get(scDim,scDim); mXXX = m_get(scDim,scDim); if (vcMean == NULL ) die("Could not malloc memory for vcMean"); if (vcShift == NULL ) die("Could not malloc memory for vcShift"); if (mSigma == NULL ) die("Could not malloc memory for mSigma"); if (mPhi == NULL ) die("Could not malloc memory for mPhi"); copyarrvec (arr,vcMean,scDim); #if DEBUG > 9 printf("\n%d scDim\n----read vcMean:%f,%f\n", scDim, vcMean->ve[0], vcMean->ve[1]); printarr(arr,scDim); printvec(vcMean); #endif sc = readfloats( Cfg[vSHIFT].data, &arr,&maxarr); if (sc != scDim ) { printf("Shift Dimension does not match process, %s!=%s\n", Cfg[vSHIFT].data, Cfg[vMEAN].data); die(""); } #if DEBUG > 9 printf("\nread shift: "); printarr(arr,scDim); #endif copyarrvec(arr,vcShift,scDim); sc = readfloats(Cfg[mSIGMA].data,&arr,&maxarr); if (sc != scDim * scDim ){ printf("%d is not enough elements for mSigma",sc); die(""); } copyarrmat(arr,mSigma,scDim,scDim); sc = readfloats(Cfg[mPHI].data,&arr,&maxarr); if (sc != scDim * scDim ){ printf("%d is not enough elements for mPhi",sc); die(""); } copyarrmat(arr,mPhi,scDim,scDim); } return 0; } /*************************************************************** Meschach Library routines, copied from meschach/{memory.c, copy.c ***************************************************************/ /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* memory.c 1.3 11/25/87 */ /* m_get -- gets an mxn matrix (in MAT form) by dynamic memory allocation */ MAT *m_get(int m, int n) { MAT *matrix; int i; if (m < 0 || n < 0) error(E_NEG,"m_get"); if ((matrix=NEW(MAT)) == (MAT *)NULL ) error(E_MEM,"m_get"); matrix->m = m; matrix->n = matrix->max_n = n; matrix->max_m = m; matrix->max_size = m*n; #ifndef SEGMENTED if ((matrix->base = NEW_A(m*n,Real)) == (Real *)NULL ) { free(matrix); error(E_MEM,"m_get"); } #else matrix->base = (Real *)NULL; #endif if ((matrix->me = (Real **)calloc((unsigned) m,sizeof(Real *))) == (Real **)NULL ) { free(matrix->base); free(matrix); error(E_MEM,"m_get"); } #ifndef SEGMENTED /* set up pointers */ for ( i=0; ime[i] = &(matrix->base[i*n]); #else for ( i = 0; i < m; i++ ) if ( (matrix->me[i]=NEW_A(n,Real)) == (Real *)NULL ) error(E_MEM,"m_get"); #endif return (matrix); } /* px_get -- gets a PERM of given 'size' by dynamic memory allocation -- Note: initialized to the identity permutation */ PERM *px_get(int size) { PERM *permute; int i; if (size < 0) error(E_NEG,"px_get"); if ((permute=NEW(PERM)) == (PERM *)NULL ) error(E_MEM,"px_get"); permute->size = permute->max_size = size; if ((permute->pe = NEW_A(size,int)) == (int *)NULL ) error(E_MEM,"px_get"); for ( i=0; ipe[i] = i; return (permute); } /* v_get -- gets a VEC of dimension 'dim' -- Note: initialized to zero */ VEC *v_get(int size) { VEC *vector; if (size < 0) error(E_NEG,"v_get"); if ((vector=NEW(VEC)) == (VEC *)NULL ) error(E_MEM,"v_get"); vector->dim = vector->max_dim = size; if ((vector->ve=NEW_A(size,Real)) == (Real *)NULL ) { free(vector); error(E_MEM,"v_get"); } return (vector); } /* m_free -- returns MAT & asoociated memory back to memory heap */ int m_free(MAT *mat) { #ifdef SEGMENTED int i; #endif if ( mat==(MAT *)NULL || (int)(mat->m) < 0 || (int)(mat->n) < 0 ) /* don't trust it */ return (-1); #ifndef SEGMENTED if ( mat->base != (Real *)NULL ) { free((char *)(mat->base)); } #else for ( i = 0; i < mat->max_m; i++ ) if ( mat->me[i] != (Real *)NULL ) { free((char *)(mat->me[i])); } #endif if ( mat->me != (Real **)NULL ) { free((char *)(mat->me)); } free((char *)mat); return (0); } /* px_free -- returns PERM & asoociated memory back to memory heap */ int px_free(px) PERM *px; { if ( px==(PERM *)NULL || (int)(px->size) < 0 ) /* don't trust it */ return (-1); if ( px->pe == (int *)NULL ) { free((char *)px); } else { free((char *)px->pe); free((char *)px); } return (0); } /* v_free -- returns VEC & asoociated memory back to memory heap */ int v_free(vec) VEC *vec; { if ( vec==(VEC *)NULL || (int)(vec->dim) < 0 ) /* don't trust it */ return (-1); if ( vec->ve == (Real *)NULL ) { free((char *)vec); } else { free((char *)vec->ve); free((char *)vec); } return (0); } /* m_resize -- returns the matrix A of size new_m x new_n; A is zeroed -- if A == NULL on entry then the effect is equivalent to m_get() */ MAT *m_resize(A,new_m,new_n) MAT *A; int new_m, new_n; { int i; int new_max_m, new_max_n, new_size, old_m, old_n; if (new_m < 0 || new_n < 0) error(E_NEG,"m_resize"); if ( ! A ) return m_get(new_m,new_n); /* nothing was changed */ if (new_m == A->m && new_n == A->n) return A; old_m = A->m; old_n = A->n; if ( new_m > A->max_m ) { /* re-allocate A->me */ A->me = RENEW(A->me,new_m,Real *); if ( ! A->me ) error(E_MEM,"m_resize"); } new_max_m = max(new_m,A->max_m); new_max_n = max(new_n,A->max_n); #ifndef SEGMENTED new_size = new_max_m*new_max_n; if ( new_size > A->max_size ) { /* re-allocate A->base */ A->base = RENEW(A->base,new_size,Real); if ( ! A->base ) error(E_MEM,"m_resize"); A->max_size = new_size; } /* now set up A->me[i] */ for ( i = 0; i < new_m; i++ ) A->me[i] = &(A->base[i*new_n]); /* now shift data in matrix */ if ( old_n > new_n ) { for ( i = 1; i < min(old_m,new_m); i++ ) MEM_COPY((char *)&(A->base[i*old_n]), (char *)&(A->base[i*new_n]), sizeof(Real)*new_n); } else if ( old_n < new_n ) { for ( i = (int)(min(old_m,new_m))-1; i > 0; i-- ) { /* copy & then zero extra space */ MEM_COPY((char *)&(A->base[i*old_n]), (char *)&(A->base[i*new_n]), sizeof(Real)*old_n); __zero__(&(A->base[i*new_n+old_n]),(new_n-old_n)); } __zero__(&(A->base[old_n]),(new_n-old_n)); A->max_n = new_n; } /* zero out the new rows.. */ for ( i = old_m; i < new_m; i++ ) __zero__(&(A->base[i*new_n]),new_n); #else if ( A->max_n < new_n ) { Real *tmp; for ( i = 0; i < A->max_m; i++ ) { if ( (tmp = RENEW(A->me[i],new_max_n,Real)) == NULL ) error(E_MEM,"m_resize"); else { A->me[i] = tmp; } } for ( i = A->max_m; i < new_max_m; i++ ) { if ( (tmp = NEW_A(new_max_n,Real)) == NULL ) error(E_MEM,"m_resize"); else { A->me[i] = tmp; } } } else if ( A->max_m < new_m ) { for ( i = A->max_m; i < new_m; i++ ) if ( (A->me[i] = NEW_A(new_max_n,Real)) == NULL ) error(E_MEM,"m_resize"); } if ( old_n < new_n ) { for ( i = 0; i < old_m; i++ ) __zero__(&(A->me[i][old_n]),new_n-old_n); } /* zero out the new rows.. */ for ( i = old_m; i < new_m; i++ ) __zero__(A->me[i],new_n); #endif A->max_m = new_max_m; A->max_n = new_max_n; A->max_size = A->max_m*A->max_n; A->m = new_m; A->n = new_n; return A; } /* px_resize -- returns the permutation px with size new_size -- px is set to the identity permutation */ PERM *px_resize(px,new_size) PERM *px; int new_size; { int i; if (new_size < 0) error(E_NEG,"px_resize"); if ( ! px ) return px_get(new_size); /* nothing is changed */ if (new_size == px->size) return px; if ( new_size > px->max_size ) { px->pe = RENEW(px->pe,new_size,int); if ( ! px->pe ) error(E_MEM,"px_resize"); px->max_size = new_size; } if ( px->size <= new_size ) /* extend permutation */ for ( i = px->size; i < new_size; i++ ) px->pe[i] = i; else for ( i = 0; i < new_size; i++ ) px->pe[i] = i; px->size = new_size; return px; } /* v_resize -- returns the vector x with dim new_dim -- x is set to the zero vector */ VEC *v_resize(x,new_dim) VEC *x; int new_dim; { if (new_dim < 0) error(E_NEG,"v_resize"); if ( ! x ) return v_get(new_dim); /* nothing is changed */ if (new_dim == x->dim) return x; if ( x->max_dim == 0 ) /* assume that it's from sub_vec */ return v_get(new_dim); if ( new_dim > x->max_dim ) { x->ve = RENEW(x->ve,new_dim,Real); if ( ! x->ve ) error(E_MEM,"v_resize"); x->max_dim = new_dim; } if ( new_dim > x->dim ) __zero__(&(x->ve[x->dim]),new_dim - x->dim); x->dim = new_dim; return x; } /* Varying number of arguments */ /* other functions of this type are in sparse.c and zmemory.c */ #ifdef ANSI_C /* To allocate memory to many arguments. The function should be called: v_get_vars(dim,&x,&y,&z,...,NULL); where int dim; VEC *x, *y, *z,...; The last argument should be NULL ! dim is the length of vectors x,y,z,... returned value is equal to the number of allocated variables Other gec_... functions are similar. */ int v_get_vars(int dim,...) { va_list ap; int i=0; VEC **par; va_start(ap, dim); while ( (par = va_arg(ap,VEC **)) ) { /* NULL ends the list*/ *par = v_get(dim); i++; } va_end(ap); return i; } int iv_get_vars(int dim,...) { va_list ap; int i=0; IVEC **par; va_start(ap, dim); while ( (par = va_arg(ap,IVEC **)) ) { /* NULL ends the list*/ *par = iv_get(dim); i++; } va_end(ap); return i; } int m_get_vars(int m,int n,...) { va_list ap; int i=0; MAT **par; va_start(ap, n); while ( (par = va_arg(ap,MAT **)) ) { /* NULL ends the list*/ *par = m_get(m,n); i++; } va_end(ap); return i; } int px_get_vars(int dim,...) { va_list ap; int i=0; PERM **par; va_start(ap, dim); while ( (par = va_arg(ap,PERM **)) ) { /* NULL ends the list*/ *par = px_get(dim); i++; } va_end(ap); return i; } /* To resize memory for many arguments. The function should be called: v_resize_vars(new_dim,&x,&y,&z,...,NULL); where int new_dim; VEC *x, *y, *z,...; The last argument should be NULL ! rdim is the resized length of vectors x,y,z,... returned value is equal to the number of allocated variables. If one of x,y,z,.. arguments is NULL then memory is allocated to this argument. Other *_resize_list() functions are similar. */ int v_resize_vars(int new_dim,...) { va_list ap; int i=0; VEC **par; va_start(ap, new_dim); while ( (par = va_arg(ap,VEC **)) ) { /* NULL ends the list*/ *par = v_resize(*par,new_dim); i++; } va_end(ap); return i; } int iv_resize_vars(int new_dim,...) { va_list ap; int i=0; IVEC **par; va_start(ap, new_dim); while ( (par = va_arg(ap,IVEC **)) ) { /* NULL ends the list*/ *par = iv_resize(*par,new_dim); i++; } va_end(ap); return i; } int m_resize_vars(int m,int n,...) { va_list ap; int i=0; MAT **par; va_start(ap, n); while ( (par = va_arg(ap,MAT **)) ) { /* NULL ends the list*/ *par = m_resize(*par,m,n); i++; } va_end(ap); return i; } int px_resize_vars(int new_dim,...) { va_list ap; int i=0; PERM **par; va_start(ap, new_dim); while ( (par = va_arg(ap,PERM **)) ) { /* NULL ends the list*/ *par = px_resize(*par,new_dim); i++; } va_end(ap); return i; } /* To deallocate memory for many arguments. The function should be called: v_free_vars(&x,&y,&z,...,NULL); where VEC *x, *y, *z,...; The last argument should be NULL ! There must be at least one not NULL argument. returned value is equal to the number of allocated variables. Returned value of x,y,z,.. is VNULL. Other *_free_list() functions are similar. */ int v_free_vars(VEC **pv,...) { va_list ap; int i=1; VEC **par; v_free(*pv); *pv = VNULL; va_start(ap, pv); while ( (par = va_arg(ap,VEC **)) ) { /* NULL ends the list*/ v_free(*par); *par = VNULL; i++; } va_end(ap); return i; } int iv_free_vars(IVEC **ipv,...) { va_list ap; int i=1; IVEC **par; iv_free(*ipv); *ipv = IVNULL; va_start(ap, ipv); while ( (par = va_arg(ap,IVEC **)) ) { /* NULL ends the list*/ iv_free(*par); *par = IVNULL; i++; } va_end(ap); return i; } int px_free_vars(PERM **vpx,...) { va_list ap; int i=1; PERM **par; px_free(*vpx); *vpx = PNULL; va_start(ap, vpx); while ( (par = va_arg(ap,PERM **)) ) { /* NULL ends the list*/ px_free(*par); *par = PNULL; i++; } va_end(ap); return i; } int m_free_vars(MAT **va,...) { va_list ap; int i=1; MAT **par; m_free(*va); *va = MNULL; va_start(ap, va); while ( (par = va_arg(ap,MAT **)) ) { /* NULL ends the list*/ m_free(*par); *par = MNULL; i++; } va_end(ap); return i; } #elif VARARGS /* old varargs is used */ /* To allocate memory to many arguments. The function should be called: v_get_vars(dim,&x,&y,&z,...,VNULL); where int dim; VEC *x, *y, *z,...; The last argument should be VNULL ! dim is the length of vectors x,y,z,... */ int v_get_vars(va_alist) va_dcl { va_list ap; int dim,i=0; VEC **par; va_start(ap); dim = va_arg(ap,int); while (par = va_arg(ap,VEC **)) { /* NULL ends the list*/ *par = v_get(dim); i++; } va_end(ap); return i; } int iv_get_vars(va_alist) va_dcl { va_list ap; int i=0, dim; IVEC **par; va_start(ap); dim = va_arg(ap,int); while (par = va_arg(ap,IVEC **)) { /* NULL ends the list*/ *par = iv_get(dim); i++; } va_end(ap); return i; } int m_get_vars(va_alist) va_dcl { va_list ap; int i=0, n, m; MAT **par; va_start(ap); m = va_arg(ap,int); n = va_arg(ap,int); while (par = va_arg(ap,MAT **)) { /* NULL ends the list*/ *par = m_get(m,n); i++; } va_end(ap); return i; } int px_get_vars(va_alist) va_dcl { va_list ap; int i=0, dim; PERM **par; va_start(ap); dim = va_arg(ap,int); while (par = va_arg(ap,PERM **)) { /* NULL ends the list*/ *par = px_get(dim); i++; } va_end(ap); return i; } /* To resize memory for many arguments. The function should be called: v_resize_vars(new_dim,&x,&y,&z,...,NULL); where int new_dim; VEC *x, *y, *z,...; The last argument should be NULL ! rdim is the resized length of vectors x,y,z,... returned value is equal to the number of allocated variables. If one of x,y,z,.. arguments is NULL then memory is allocated to this argument. Other *_resize_list() functions are similar. */ int v_resize_vars(va_alist) va_dcl { va_list ap; int i=0, new_dim; VEC **par; va_start(ap); new_dim = va_arg(ap,int); while (par = va_arg(ap,VEC **)) { /* NULL ends the list*/ *par = v_resize(*par,new_dim); i++; } va_end(ap); return i; } int iv_resize_vars(va_alist) va_dcl { va_list ap; int i=0, new_dim; IVEC **par; va_start(ap); new_dim = va_arg(ap,int); while (par = va_arg(ap,IVEC **)) { /* NULL ends the list*/ *par = iv_resize(*par,new_dim); i++; } va_end(ap); return i; } int m_resize_vars(va_alist) va_dcl { va_list ap; int i=0, m, n; MAT **par; va_start(ap); m = va_arg(ap,int); n = va_arg(ap,int); while (par = va_arg(ap,MAT **)) { /* NULL ends the list*/ *par = m_resize(*par,m,n); i++; } va_end(ap); return i; } int px_resize_vars(va_alist) va_dcl { va_list ap; int i=0, new_dim; PERM **par; va_start(ap); new_dim = va_arg(ap,int); while (par = va_arg(ap,PERM **)) { /* NULL ends the list*/ *par = px_resize(*par,new_dim); i++; } va_end(ap); return i; } /* To deallocate memory for many arguments. The function should be called: v_free_vars(&x,&y,&z,...,NULL); where VEC *x, *y, *z,...; The last argument should be NULL ! returned value is equal to the number of allocated variables. Returned value of x,y,z,.. is VNULL. Other *_free_list() functions are similar. */ int v_free_vars(va_alist) va_dcl { va_list ap; int i=0; VEC **par; va_start(ap); while (par = va_arg(ap,VEC **)) { /* NULL ends the list*/ v_free(*par); *par = VNULL; i++; } va_end(ap); return i; } int iv_free_vars(va_alist) va_dcl { va_list ap; int i=0; IVEC **par; va_start(ap); while (par = va_arg(ap,IVEC **)) { /* NULL ends the list*/ iv_free(*par); *par = IVNULL; i++; } va_end(ap); return i; } int px_free_vars(va_alist) va_dcl { va_list ap; int i=0; PERM **par; va_start(ap); while (par = va_arg(ap,PERM **)) { /* NULL ends the list*/ px_free(*par); *par = PNULL; i++; } va_end(ap); return i; } int m_free_vars(va_alist) va_dcl { va_list ap; int i=0; MAT **par; va_start(ap); while (par = va_arg(ap,MAT **)) { /* NULL ends the list*/ m_free(*par); *par = MNULL; i++; } va_end(ap); return i; } #endif /* VARARGS */ /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* File containing routines for symmetric eigenvalue problems */ #define SQRT2 1.4142135623730949 #define sgn(x) ( (x) >= 0 ? 1 : -1 ) /* trieig -- finds eigenvalues of symmetric tridiagonal matrices -- matrix represented by a pair of vectors a (diag entries) and b (sub- & super-diag entries) -- eigenvalues in a on return */ VEC *trieig(VEC* a,VEC* b,MAT*Q) { int i, i_min, i_max, n, split; Real *a_ve, *b_ve; Real b_sqr, bk, ak1, bk1, ak2, bk2, z; Real c, c2, cs, s, s2, d, mu; if ( ! a || ! b ) error(E_NULL,"trieig"); if ( a->dim != b->dim + 1 || ( Q && Q->m != a->dim ) ) error(E_SIZES,"trieig"); if ( Q && Q->m != Q->n ) error(E_SQUARE,"trieig"); n = a->dim; a_ve = a->ve; b_ve = b->ve; i_min = 0; while ( i_min < n ) /* outer while loop */ { /* find i_max to suit; submatrix i_min..i_max should be irreducible */ i_max = n-1; for ( i = i_min; i < n-1; i++ ) if ( b_ve[i] == 0.0 ) { i_max = i; break; } if ( i_max <= i_min ) { /* printf("# i_min = %d, i_max = %d\n",i_min,i_max); */ i_min = i_max + 1; continue; /* outer while loop */ } /* printf("# i_min = %d, i_max = %d\n",i_min,i_max); */ /* repeatedly perform QR method until matrix splits */ split = FALSE; while ( ! split ) /* inner while loop */ { /* find Wilkinson shift */ d = (a_ve[i_max-1] - a_ve[i_max])/2; b_sqr = b_ve[i_max-1]*b_ve[i_max-1]; mu = a_ve[i_max] - b_sqr/(d + sgn(d)*sqrt(d*d+b_sqr)); /* printf("# Wilkinson shift = %g\n",mu); */ /* initial Givens' rotation */ givens(a_ve[i_min]-mu,b_ve[i_min],&c,&s); s = -s; /* printf("# c = %g, s = %g\n",c,s); */ if ( fabs(c) < SQRT2 ) { c2 = c*c; s2 = 1-c2; } else { s2 = s*s; c2 = 1-s2; } cs = c*s; ak1 = c2*a_ve[i_min]+s2*a_ve[i_min+1]-2*cs*b_ve[i_min]; bk1 = cs*(a_ve[i_min]-a_ve[i_min+1]) + (c2-s2)*b_ve[i_min]; ak2 = s2*a_ve[i_min]+c2*a_ve[i_min+1]+2*cs*b_ve[i_min]; bk2 = ( i_min < i_max-1 ) ? c*b_ve[i_min+1] : 0.0; z = ( i_min < i_max-1 ) ? -s*b_ve[i_min+1] : 0.0; a_ve[i_min] = ak1; a_ve[i_min+1] = ak2; b_ve[i_min] = bk1; if ( i_min < i_max-1 ) b_ve[i_min+1] = bk2; if ( Q ) rot_cols(Q,i_min,i_min+1,c,-s,Q); /* printf("# z = %g\n",z); */ /* printf("# a [temp1] =\n"); v_output(a); */ /* printf("# b [temp1] =\n"); v_output(b); */ for ( i = i_min+1; i < i_max; i++ ) { /* get Givens' rotation for sub-block -- k == i-1 */ givens(b_ve[i-1],z,&c,&s); s = -s; /* printf("# c = %g, s = %g\n",c,s); */ /* perform Givens' rotation on sub-block */ if ( fabs(c) < SQRT2 ) { c2 = c*c; s2 = 1-c2; } else { s2 = s*s; c2 = 1-s2; } cs = c*s; bk = c*b_ve[i-1] - s*z; ak1 = c2*a_ve[i]+s2*a_ve[i+1]-2*cs*b_ve[i]; bk1 = cs*(a_ve[i]-a_ve[i+1]) + (c2-s2)*b_ve[i]; ak2 = s2*a_ve[i]+c2*a_ve[i+1]+2*cs*b_ve[i]; bk2 = ( i+1 < i_max ) ? c*b_ve[i+1] : 0.0; z = ( i+1 < i_max ) ? -s*b_ve[i+1] : 0.0; a_ve[i] = ak1; a_ve[i+1] = ak2; b_ve[i] = bk1; if ( i < i_max-1 ) b_ve[i+1] = bk2; if ( i > i_min ) b_ve[i-1] = bk; if ( Q ) rot_cols(Q,i,i+1,c,-s,Q); /* printf("# a [temp2] =\n"); v_output(a); */ /* printf("# b [temp2] =\n"); v_output(b); */ } /* test to see if matrix should be split */ for ( i = i_min; i < i_max; i++ ) if ( fabs(b_ve[i]) < MACHEPS* (fabs(a_ve[i])+fabs(a_ve[i+1])) ) { b_ve[i] = 0.0; split = TRUE; } /* printf("# a =\n"); v_output(a); */ /* printf("# b =\n"); v_output(b); */ } } return a; } /* symmeig -- computes eigenvalues of a dense symmetric matrix -- A **must** be symmetric on entry -- eigenvalues stored in out -- Q contains orthogonal matrix of eigenvectors -- returns vector of eigenvalues */ VEC *symmeig(MAT* A,MAT* Q,VEC* out) { int i; static MAT *tmp = MNULL; static VEC *b = VNULL, *diag = VNULL, *beta = VNULL; if ( ! A ) error(E_NULL,"symmeig"); if ( A->m != A->n ) error(E_SQUARE,"symmeig"); if ( ! out || out->dim != A->m ) out = v_resize(out,A->m); tmp = m_copy(A,tmp); b = v_resize(b,A->m - 1); diag = v_resize(diag,(int)A->m); beta = v_resize(beta,(int)A->m); MEM_STAT_REG(tmp,TYPE_MAT); MEM_STAT_REG(b,TYPE_VEC); MEM_STAT_REG(diag,TYPE_VEC); MEM_STAT_REG(beta,TYPE_VEC); Hfactor(tmp,diag,beta); if ( Q ) makeHQ(tmp,diag,beta,Q); for ( i = 0; i < A->m - 1; i++ ) { out->ve[i] = tmp->me[i][i]; b->ve[i] = tmp->me[i][i+1]; } out->ve[i] = tmp->me[i][i]; trieig(out,b,Q); return out; } /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* matop.c 1.3 11/25/87 */ /* m_add -- matrix addition -- may be in-situ */ MAT *m_add(mat1,mat2,out) MAT *mat1,*mat2,*out; { int m,n,i; if ( mat1==(MAT *)NULL || mat2==(MAT *)NULL ) error(E_NULL,"m_add"); if ( mat1->m != mat2->m || mat1->n != mat2->n ) error(E_SIZES,"m_add"); if ( out==(MAT *)NULL || out->m != mat1->m || out->n != mat1->n ) out = m_resize(out,mat1->m,mat1->n); m = mat1->m; n = mat1->n; for ( i=0; ime[i],mat2->me[i],out->me[i],n); /************************************************** for ( j=0; jme[i][j] = mat1->me[i][j]+mat2->me[i][j]; **************************************************/ } return (out); } /* m_sub -- matrix subtraction -- may be in-situ */ MAT *m_sub(mat1,mat2,out) MAT *mat1,*mat2,*out; { int m,n,i; if ( mat1==(MAT *)NULL || mat2==(MAT *)NULL ) error(E_NULL,"m_sub"); if ( mat1->m != mat2->m || mat1->n != mat2->n ) error(E_SIZES,"m_sub"); if ( out==(MAT *)NULL || out->m != mat1->m || out->n != mat1->n ) out = m_resize(out,mat1->m,mat1->n); m = mat1->m; n = mat1->n; for ( i=0; ime[i],mat2->me[i],out->me[i],n); /************************************************** for ( j=0; jme[i][j] = mat1->me[i][j]-mat2->me[i][j]; **************************************************/ } return (out); } /* m_mlt -- matrix-matrix multiplication */ MAT *m_mlt(A,B,OUT) MAT *A,*B,*OUT; { int i, /* j, */ k, m, n, p; Real **A_v, **B_v /*, *B_row, *OUT_row, sum, tmp */; if ( A==(MAT *)NULL || B==(MAT *)NULL ) error(E_NULL,"m_mlt"); if ( A->n != B->m ) error(E_SIZES,"m_mlt"); if ( A == OUT || B == OUT ) error(E_INSITU,"m_mlt"); m = A->m; n = A->n; p = B->n; A_v = A->me; B_v = B->me; if ( OUT==(MAT *)NULL || OUT->m != A->m || OUT->n != B->n ) OUT = m_resize(OUT,A->m,B->n); /**************************************************************** for ( i=0; ime[i][j] = sum; } ****************************************************************/ m_zero(OUT); for ( i=0; ime[i],B_v[k],A_v[i][k],p); /************************************************** B_row = B_v[k]; OUT_row = OUT->me[i]; for ( j=0; jn != B->n ) error(E_SIZES,"mmtr_mlt"); if ( ! OUT || OUT->m != A->m || OUT->n != B->m ) OUT = m_resize(OUT,A->m,B->m); limit = A->n; for ( i = 0; i < A->m; i++ ) for ( j = 0; j < B->m; j++ ) { OUT->me[i][j] = __ip__(A->me[i],B->me[j],(int)limit); /************************************************** sum = 0.0; A_row = A->me[i]; B_row = B->me[j]; for ( k = 0; k < limit; k++ ) sum += (*A_row++)*(*B_row++); OUT->me[i][j] = sum; **************************************************/ } return OUT; } /* mtrm_mlt -- matrix transposed-matrix multiplication -- A^T.B is returned, result stored in OUT */ MAT *mtrm_mlt(A,B,OUT) MAT *A, *B, *OUT; { int i, k, limit; /* Real *B_row, *OUT_row, multiplier; */ if ( ! A || ! B ) error(E_NULL,"mmtr_mlt"); if ( A == OUT || B == OUT ) error(E_INSITU,"mtrm_mlt"); if ( A->m != B->m ) error(E_SIZES,"mmtr_mlt"); if ( ! OUT || OUT->m != A->n || OUT->n != B->n ) OUT = m_resize(OUT,A->n,B->n); limit = B->n; m_zero(OUT); for ( k = 0; k < A->m; k++ ) for ( i = 0; i < A->n; i++ ) { if ( A->me[k][i] != 0.0 ) __mltadd__(OUT->me[i],B->me[k],A->me[k][i],limit); /************************************************** multiplier = A->me[k][i]; OUT_row = OUT->me[i]; B_row = B->me[k]; for ( j = 0; j < limit; j++ ) *(OUT_row++) += multiplier*(*B_row++); **************************************************/ } return OUT; } /* mv_mlt -- matrix-vector multiplication -- Note: b is treated as a column vector */ VEC *mv_mlt(A,b,out) MAT *A; VEC *b,*out; { int i, m, n; Real **A_v, *b_v /*, *A_row */; /* register Real sum; */ if ( A==(MAT *)NULL || b==(VEC *)NULL ) error(E_NULL,"mv_mlt"); if ( A->n != b->dim ) error(E_SIZES,"mv_mlt"); if ( b == out ) error(E_INSITU,"mv_mlt"); if ( out == (VEC *)NULL || out->dim != A->m ) out = v_resize(out,A->m); m = A->m; n = A->n; A_v = A->me; b_v = b->ve; for ( i=0; ive[i] = __ip__(A_v[i],b_v,(int)n); /************************************************** A_row = A_v[i]; b_v = b->ve; for ( j=0; jve[i] = sum; **************************************************/ } return out; } /* sm_mlt -- scalar-matrix multiply -- may be in-situ */ MAT *sm_mlt(scalar,matrix,out) double scalar; MAT *matrix,*out; { int m,n,i; if ( matrix==(MAT *)NULL ) error(E_NULL,"sm_mlt"); if ( out==(MAT *)NULL || out->m != matrix->m || out->n != matrix->n ) out = m_resize(out,matrix->m,matrix->n); m = matrix->m; n = matrix->n; for ( i=0; ime[i],(double)scalar,out->me[i],n); /************************************************** for ( j=0; jme[i][j] = scalar*matrix->me[i][j]; **************************************************/ return (out); } /* vm_mlt -- vector-matrix multiplication -- Note: b is treated as a row vector */ VEC *vm_mlt(A,b,out) MAT *A; VEC *b,*out; { int j,m,n; /* Real sum,**A_v,*b_v; */ if ( A==(MAT *)NULL || b==(VEC *)NULL ) error(E_NULL,"vm_mlt"); if ( A->m != b->dim ) error(E_SIZES,"vm_mlt"); if ( b == out ) error(E_INSITU,"vm_mlt"); if ( out == (VEC *)NULL || out->dim != A->n ) out = v_resize(out,A->n); m = A->m; n = A->n; v_zero(out); for ( j = 0; j < m; j++ ) if ( b->ve[j] != 0.0 ) __mltadd__(out->ve,A->me[j],b->ve[j],n); /************************************************** A_v = A->me; b_v = b->ve; for ( j=0; jve[j] = sum; } **************************************************/ return out; } /* m_transp -- transpose matrix */ MAT *m_transp(in,out) MAT *in, *out; { int i, j; int in_situ; Real tmp; if ( in == (MAT *)NULL ) error(E_NULL,"m_transp"); if ( in == out && in->n != in->m ) error(E_INSITU2,"m_transp"); in_situ = ( in == out ); if ( out == (MAT *)NULL || out->m != in->n || out->n != in->m ) out = m_resize(out,in->n,in->m); if ( ! in_situ ) for ( i = 0; i < in->m; i++ ) for ( j = 0; j < in->n; j++ ) out->me[j][i] = in->me[i][j]; else for ( i = 1; i < in->m; i++ ) for ( j = 0; j < i; j++ ) { tmp = in->me[i][j]; in->me[i][j] = in->me[j][i]; in->me[j][i] = tmp; } return out; } /* swap_rows -- swaps rows i and j of matrix A upto column lim */ MAT *swap_rows(A,i,j,lo,hi) MAT *A; int i, j, lo, hi; { int k; Real **A_me, tmp; if ( ! A ) error(E_NULL,"swap_rows"); if ( i < 0 || j < 0 || i >= A->m || j >= A->m ) error(E_SIZES,"swap_rows"); lo = max(0,lo); hi = min(hi,A->n-1); A_me = A->me; for ( k = lo; k <= hi; k++ ) { tmp = A_me[k][i]; A_me[k][i] = A_me[k][j]; A_me[k][j] = tmp; } return A; } /* swap_cols -- swap columns i and j of matrix A upto row lim */ MAT *swap_cols(A,i,j,lo,hi) MAT *A; int i, j, lo, hi; { int k; Real **A_me, tmp; if ( ! A ) error(E_NULL,"swap_cols"); if ( i < 0 || j < 0 || i >= A->n || j >= A->n ) error(E_SIZES,"swap_cols"); lo = max(0,lo); hi = min(hi,A->m-1); A_me = A->me; for ( k = lo; k <= hi; k++ ) { tmp = A_me[i][k]; A_me[i][k] = A_me[j][k]; A_me[j][k] = tmp; } return A; } /* ms_mltadd -- matrix-scalar multiply and add -- may be in situ -- returns out == A1 + s*A2 */ MAT *ms_mltadd(A1,A2,s,out) MAT *A1, *A2, *out; double s; { /* register Real *A1_e, *A2_e, *out_e; */ /* register int j; */ int i, m, n; if ( ! A1 || ! A2 ) error(E_NULL,"ms_mltadd"); if ( A1->m != A2->m || A1->n != A2->n ) error(E_SIZES,"ms_mltadd"); if ( s == 0.0 ) return m_copy(A1,out); if ( s == 1.0 ) return m_add(A1,A2,out); tracecatch(out = m_copy(A1,out),"ms_mltadd"); m = A1->m; n = A1->n; for ( i = 0; i < m; i++ ) { __mltadd__(out->me[i],A2->me[i],s,n); /************************************************** A1_e = A1->me[i]; A2_e = A2->me[i]; out_e = out->me[i]; for ( j = 0; j < n; j++ ) out_e[j] = A1_e[j] + s*A2_e[j]; **************************************************/ } return out; } /* mv_mltadd -- matrix-vector multiply and add -- may not be in situ -- returns out == v1 + alpha*A*v2 */ VEC *mv_mltadd(v1,v2,A,alpha,out) VEC *v1, *v2, *out; MAT *A; double alpha; { /* register int j; */ int i, m, n; Real *v2_ve, *out_ve; if ( ! v1 || ! v2 || ! A ) error(E_NULL,"mv_mltadd"); if ( out == v2 ) error(E_INSITU,"mv_mltadd"); if ( v1->dim != A->m || v2->dim != A-> n ) error(E_SIZES,"mv_mltadd"); tracecatch(out = v_copy(v1,out),"mv_mltadd"); v2_ve = v2->ve; out_ve = out->ve; m = A->m; n = A->n; if ( alpha == 0.0 ) return out; for ( i = 0; i < m; i++ ) { out_ve[i] += alpha*__ip__(A->me[i],v2_ve,(int)n); /************************************************** A_e = A->me[i]; sum = 0.0; for ( j = 0; j < n; j++ ) sum += A_e[j]*v2_ve[j]; out_ve[i] = v1->ve[i] + alpha*sum; **************************************************/ } return out; } /* vm_mltadd -- vector-matrix multiply and add -- may not be in situ -- returns out' == v1' + v2'*A */ VEC *vm_mltadd(v1,v2,A,alpha,out) VEC *v1, *v2, *out; MAT *A; double alpha; { int /* i, */ j, m, n; Real tmp, /* *A_e, */ *out_ve; if ( ! v1 || ! v2 || ! A ) error(E_NULL,"vm_mltadd"); if ( v2 == out ) error(E_INSITU,"vm_mltadd"); if ( v1->dim != A->n || A->m != v2->dim ) error(E_SIZES,"vm_mltadd"); tracecatch(out = v_copy(v1,out),"vm_mltadd"); out_ve = out->ve; m = A->m; n = A->n; for ( j = 0; j < m; j++ ) { tmp = v2->ve[j]*alpha; if ( tmp != 0.0 ) __mltadd__(out_ve,A->me[j],tmp,n); /************************************************** A_e = A->me[j]; for ( i = 0; i < n; i++ ) out_ve[i] += A_e[i]*tmp; **************************************************/ } return out; } /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* ivecop.c */ static char line[MAXLINE]; /* iv_get -- get integer vector -- see also memory.c */ IVEC *iv_get(dim) int dim; { IVEC *iv; /* int i; */ if (dim < 0) error(E_NEG,"iv_get"); if ((iv=NEW(IVEC)) == IVNULL ) error(E_MEM,"iv_get"); iv->dim = iv->max_dim = dim; if ((iv->ive = NEW_A(dim,int)) == (int *)NULL ) error(E_MEM,"iv_get"); return (iv); } /* iv_free -- returns iv & asoociated memory back to memory heap */ int iv_free(iv) IVEC *iv; { if ( iv==IVNULL || iv->dim > MAXDIM ) /* don't trust it */ return (-1); if ( iv->ive == (int *)NULL ) { free((char *)iv); } else { free((char *)iv->ive); free((char *)iv); } return (0); } /* iv_resize -- returns the IVEC with dimension new_dim -- iv is set to the zero vector */ IVEC *iv_resize(iv,new_dim) IVEC *iv; int new_dim; { int i; if (new_dim < 0) error(E_NEG,"iv_resize"); if ( ! iv ) return iv_get(new_dim); if (new_dim == iv->dim) return iv; if ( new_dim > iv->max_dim ) { iv->ive = RENEW(iv->ive,new_dim,int); if ( ! iv->ive ) error(E_MEM,"iv_resize"); iv->max_dim = new_dim; } if ( iv->dim <= new_dim ) for ( i = iv->dim; i < new_dim; i++ ) iv->ive[i] = 0; iv->dim = new_dim; return iv; } /* iv_copy -- copy integer vector in to out -- out created/resized if necessary */ IVEC *iv_copy(in,out) IVEC *in, *out; { int i; if ( ! in ) error(E_NULL,"iv_copy"); out = iv_resize(out,in->dim); for ( i = 0; i < in->dim; i++ ) out->ive[i] = in->ive[i]; return out; } /* iv_move -- move selected pieces of an IVEC -- moves the length dim0 subvector with initial index i0 to the corresponding subvector of out with initial index i1 -- out is resized if necessary */ IVEC *iv_move(in,i0,dim0,out,i1) IVEC *in, *out; int i0, dim0, i1; { if ( ! in ) error(E_NULL,"iv_move"); if ( i0 < 0 || dim0 < 0 || i1 < 0 || i0+dim0 > in->dim ) error(E_BOUNDS,"iv_move"); if ( (! out) || i1+dim0 > out->dim ) out = iv_resize(out,i1+dim0); MEM_COPY(&(in->ive[i0]),&(out->ive[i1]),dim0*sizeof(int)); return out; } /* iv_add -- integer vector addition -- may be in-situ */ IVEC *iv_add(iv1,iv2,out) IVEC *iv1,*iv2,*out; { int i; int *out_ive, *iv1_ive, *iv2_ive; if ( iv1==IVNULL || iv2==IVNULL ) error(E_NULL,"iv_add"); if ( iv1->dim != iv2->dim ) error(E_SIZES,"iv_add"); if ( out==IVNULL || out->dim != iv1->dim ) out = iv_resize(out,iv1->dim); out_ive = out->ive; iv1_ive = iv1->ive; iv2_ive = iv2->ive; for ( i = 0; i < iv1->dim; i++ ) out_ive[i] = iv1_ive[i] + iv2_ive[i]; return (out); } /* iv_sub -- integer vector addition -- may be in-situ */ IVEC *iv_sub(iv1,iv2,out) IVEC *iv1,*iv2,*out; { int i; int *out_ive, *iv1_ive, *iv2_ive; if ( iv1==IVNULL || iv2==IVNULL ) error(E_NULL,"iv_sub"); if ( iv1->dim != iv2->dim ) error(E_SIZES,"iv_sub"); if ( out==IVNULL || out->dim != iv1->dim ) out = iv_resize(out,iv1->dim); out_ive = out->ive; iv1_ive = iv1->ive; iv2_ive = iv2->ive; for ( i = 0; i < iv1->dim; i++ ) out_ive[i] = iv1_ive[i] - iv2_ive[i]; return (out); } /* iv_foutput -- print a representation of iv on stream fp */ void iv_foutput(fp,iv) FILE *fp; IVEC *iv; { int i; fprintf(fp,"IntVector: "); if ( iv == IVNULL ) { fprintf(fp,"**** NULL ****\n"); return; } fprintf(fp,"dim: %d\n",iv->dim); for ( i = 0; i < iv->dim; i++ ) { if ( (i+1) % 8 ) fprintf(fp,"%8d ",iv->ive[i]); else fprintf(fp,"%8d\n",iv->ive[i]); } if ( i % 8 ) fprintf(fp,"\n"); } /* iv_finput -- input integer vector from stream fp */ IVEC *iv_finput(FILE* fp, IVEC* x) { IVEC *iiv_finput(FILE*,IVEC*),*biv_finput(FILE*,IVEC*); if ( isatty((int)fileno(fp)) ) return iiv_finput(fp,x); else return biv_finput(fp,x); } /* iiv_finput -- interactive input of IVEC iv */ IVEC *iiv_finput(FILE* fp, IVEC* iv) { int i,dim,dynamic; /* dynamic set if memory allocated here */ /* get dimension */ if ( iv != (IVEC *)NULL && iv->dimdim; dynamic = FALSE; } else { dynamic = TRUE; do { fprintf(stderr,"IntVector: dim: "); if ( fgets(line,MAXLINE,fp)==NULL ) error(E_INPUT,"iiv_finput"); } while ( sscanf(line,"%u",&dim)<1 || dim>MAXDIM ); iv = iv_get(dim); } /* input elements */ for ( i=0; iive[i]); if ( fgets(line,MAXLINE,fp)==NULL ) error(E_INPUT,"iiv_finput"); if ( (*line == 'b' || *line == 'B') && i > 0 ) { i--; dynamic = FALSE; goto redo; } if ( (*line == 'f' || *line == 'F') && i < dim-1 ) { i++; dynamic = FALSE; goto redo; } } while ( *line=='\0' || sscanf(line,"%d",&iv->ive[i]) < 1 ); return (iv); } /* biv_finput -- batch-file input of IVEC iv */ IVEC *biv_finput(FILE* fp, IVEC* iv) { int i,dim; int io_code; /* get dimension */ skipjunk(fp); if ((io_code=fscanf(fp," IntVector: dim:%u",&dim)) < 1 || dim>MAXDIM ) error(io_code==EOF ? 7 : 6,"biv_finput"); /* allocate memory if necessary */ if ( iv==(IVEC *)NULL || iv->dimive[i])) < 1 ) error(io_code==EOF ? 7 : 6,"biv_finput"); return (iv); } /* iv_dump -- dumps all the contents of IVEC iv onto stream fp */ void iv_dump(fp,iv) FILE*fp; IVEC*iv; { int i; fprintf(fp,"IntVector: "); if ( ! iv ) { fprintf(fp,"**** NULL ****\n"); return; } fprintf(fp,"dim: %d, max_dim: %d\n",iv->dim,iv->max_dim); fprintf(fp,"ive @ 0x%lx\n",(long)(iv->ive)); for ( i = 0; i < iv->max_dim; i++ ) { if ( (i+1) % 8 ) fprintf(fp,"%8d ",iv->ive[i]); else fprintf(fp,"%8d\n",iv->ive[i]); } if ( i % 8 ) fprintf(fp,"\n"); } #define MAX_STACK 60 /* iv_sort -- sorts vector x, and generates permutation that gives the order of the components; x = [1.3, 3.7, 0.5] -> [0.5, 1.3, 3.7] and the permutation is order = [2, 0, 1]. -- if order is NULL on entry then it is ignored -- the sorted vector x is returned */ IVEC *iv_sort(x, order) IVEC *x; PERM *order; { int *x_ive, tmp, v; /* int *order_pe; */ int dim, i, j, l, r, tmp_i; int stack[MAX_STACK], sp; if ( ! x ) error(E_NULL,"v_sort"); if ( order != PNULL && order->size != x->dim ) order = px_resize(order, x->dim); x_ive = x->ive; dim = x->dim; if ( order != PNULL ) px_ident(order); if ( dim <= 1 ) return x; /* using quicksort algorithm in Sedgewick, "Algorithms in C", Ch. 9, pp. 118--122 (1990) */ sp = 0; l = 0; r = dim-1; v = x_ive[0]; for ( ; ; ) { while ( r > l ) { /* "i = partition(x_ive,l,r);" */ v = x_ive[r]; i = l-1; j = r; for ( ; ; ) { while ( x_ive[++i] < v ) ; while ( x_ive[--j] > v ) ; if ( i >= j ) break; tmp = x_ive[i]; x_ive[i] = x_ive[j]; x_ive[j] = tmp; if ( order != PNULL ) { tmp_i = order->pe[i]; order->pe[i] = order->pe[j]; order->pe[j] = tmp_i; } } tmp = x_ive[i]; x_ive[i] = x_ive[r]; x_ive[r] = tmp; if ( order != PNULL ) { tmp_i = order->pe[i]; order->pe[i] = order->pe[r]; order->pe[r] = tmp_i; } if ( i-l > r-i ) { stack[sp++] = l; stack[sp++] = i-1; l = i+1; } else { stack[sp++] = i+1; stack[sp++] = r; r = i-1; } } /* recursion elimination */ if ( sp == 0 ) break; r = stack[--sp]; l = stack[--sp]; } return x; } /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* Files for matrix computations Givens operations file. Contains routines for calculating and applying givens rotations for/to vectors and also to matrices by row and by column. */ /* givens.c 1.2 11/25/87 */ /* givens -- returns c,s parameters for Givens rotation to eliminate y in the vector [ x y ]' */ void givens(x,y,c,s) double x,y; Real *c,*s; { Real norm; norm = sqrt(x*x+y*y); if ( norm == 0.0 ) { *c = 1.0; *s = 0.0; } /* identity */ else { *c = x/norm; *s = y/norm; } } /* rot_vec -- apply Givens rotation to x's i & k components */ VEC *rot_vec(x,i,k,c,s,out) VEC *x,*out; int i,k; double c,s; { Real temp; if ( x==VNULL ) error(E_NULL,"rot_vec"); if ( i >= x->dim || k >= x->dim ) error(E_RANGE,"rot_vec"); out = v_copy(x,out); /* temp = c*out->ve[i] + s*out->ve[k]; */ temp = c*v_entry(out,i) + s*v_entry(out,k); /* out->ve[k] = -s*out->ve[i] + c*out->ve[k]; */ v_set_val(out,k,-s*v_entry(out,i)+c*v_entry(out,k)); /* out->ve[i] = temp; */ v_set_val(out,i,temp); return (out); } /* rot_rows -- premultiply mat by givens rotation described by c,s */ MAT *rot_rows(mat,i,k,c,s,out) MAT *mat,*out; int i,k; double c,s; { int j; Real temp; if ( mat==(MAT *)NULL ) error(E_NULL,"rot_rows"); if ( i >= mat->m || k >= mat->m ) error(E_RANGE,"rot_rows"); out = m_copy(mat,out); for ( j=0; jn; j++ ) { /* temp = c*out->me[i][j] + s*out->me[k][j]; */ temp = c*m_entry(out,i,j) + s*m_entry(out,k,j); /* out->me[k][j] = -s*out->me[i][j] + c*out->me[k][j]; */ m_set_val(out,k,j, -s*m_entry(out,i,j) + c*m_entry(out,k,j)); /* out->me[i][j] = temp; */ m_set_val(out,i,j, temp); } return (out); } /* rot_cols -- postmultiply mat by givens rotation described by c,s */ MAT *rot_cols(mat,i,k,c,s,out) MAT *mat,*out; int i,k; double c,s; { int j; Real temp; if ( mat==(MAT *)NULL ) error(E_NULL,"rot_cols"); if ( i < 0 || i >= mat->n || k < 0 || k >= mat->n ) error(E_RANGE,"rot_cols"); out = m_copy(mat,out); for ( j=0; jm; j++ ) { /* temp = c*out->me[j][i] + s*out->me[j][k]; */ temp = c*m_entry(out,j,i) + s*m_entry(out,j,k); /* out->me[j][k] = -s*out->me[j][i] + c*out->me[j][k]; */ m_set_val(out,j,k, -s*m_entry(out,j,i) + c*m_entry(out,j,k)); /* out->me[j][i] = temp; */ m_set_val(out,j,i,temp); } return (out); } /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* mem_stat.c 6/09/93 */ /* Deallocation of static arrays */ /* global variable */ extern MEM_CONNECT mem_connect[MEM_CONNECT_MAX_LISTS]; /* names of types */ static char *mem_type_names[] = { "MAT", "BAND", "PERM", "VEC", "IVEC" #ifdef SPARSE ,"ITER", "SPROW", "SPMAT" #endif #ifdef COMPLEX ,"ZVEC", "ZMAT" #endif }; #define MEM_NUM_STD_TYPES (sizeof(mem_type_names)/sizeof(mem_type_names[0])) /* local array for keeping track of memory */ static MEM_ARRAY mem_info_sum[MEM_NUM_STD_TYPES]; /* for freeing various types */ static int (*mem_free_funcs[MEM_NUM_STD_TYPES])(void*) = { (int(*)(void*)) m_free, (int(*)(void*)) bd_free, (int(*)(void*)) px_free, (int(*)(void*)) v_free, (int(*)(void*)) iv_free #ifdef SPARSE ,iter_free, sprow_free, sp_free #endif #ifdef COMPLEX ,zv_free, zm_free #endif }; /* it is a global variable for passing pointers to local arrays defined here */ MEM_CONNECT mem_connect[MEM_CONNECT_MAX_LISTS] = { { mem_type_names, mem_free_funcs, MEM_NUM_STD_TYPES, mem_info_sum } }; /* local type */ typedef struct { void **var; /* for &A, where A is a pointer */ int type; /* type of A */ int mark; /* what mark is chosen */ } MEM_STAT_STRUCT; /* local variables */ /* how many marks are used */ static int mem_stat_mark_many = 0; /* current mark */ static int mem_stat_mark_curr = 0; static MEM_STAT_STRUCT mem_stat_var[MEM_HASHSIZE]; /* array of indices (+1) to mem_stat_var */ static unsigned int mem_hash_idx[MEM_HASHSIZE]; /* points to the first unused element in mem_hash_idx */ static unsigned int mem_hash_idx_end = 0; /* hashing function */ static unsigned int mem_hash(void ** ptr) { unsigned long lp = (unsigned long)ptr; return (lp % MEM_HASHSIZE); } /* look for a place in mem_stat_var */ static int mem_lookup(void **var) { int k, j; k = mem_hash(var); if (mem_stat_var[k].var == var) { return -1; } else if (mem_stat_var[k].var == NULL) { return k; } else { /* look for an empty place */ j = k; while (mem_stat_var[j].var != var && j < MEM_HASHSIZE && mem_stat_var[j].var != NULL) j++; if (mem_stat_var[j].var == NULL) return j; else if (mem_stat_var[j].var == var) return -1; else { /* if (j == MEM_HASHSIZE) */ j = 0; while (mem_stat_var[j].var != var && j < k && mem_stat_var[j].var != NULL) j++; if (mem_stat_var[j].var == NULL) return j; else if (mem_stat_var[j].var == var) return -1; else { /* if (j == k) */ fprintf(stderr, "\n WARNING !!! static memory: mem_stat_var is too small\n"); fprintf(stderr, " Increase MEM_HASHSIZE in file: %s (currently = %d)\n\n", MEM_HASHSIZE_FILE, MEM_HASHSIZE); if ( !isatty(fileno(stdout)) ) { fprintf(stdout, "\n WARNING !!! static memory: mem_stat_var is too small\n"); fprintf(stdout, " Increase MEM_HASHSIZE in file: %s (currently = %d)\n\n", MEM_HASHSIZE_FILE, MEM_HASHSIZE); } error(E_MEM,"mem_lookup"); } } } return -1; } /* register static variables; Input arguments: var - variable to be registered, type - type of this variable; list - list of types returned value < 0 --> error, returned value == 0 --> not registered, returned value >= 0 --> registered with this mark; */ int mem_stat_reg_list(var,type,list) void **var; int type,list; { int n; if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS ) return -1; if (mem_stat_mark_curr == 0) return 0; /* not registered */ if (var == NULL) return -1; /* error */ if ( type < 0 || (unsigned) type >= mem_connect[list].ntypes || mem_connect[list].free_funcs[type] == NULL ) { warning(WARN_WRONG_TYPE,"mem_stat_reg_list"); return -1; } if ((n = mem_lookup(var)) >= 0) { mem_stat_var[n].var = var; mem_stat_var[n].mark = mem_stat_mark_curr; mem_stat_var[n].type = type; /* save n+1, not n */ mem_hash_idx[mem_hash_idx_end++] = n+1; } return mem_stat_mark_curr; } /* set a mark; Input argument: mark - positive number denoting a mark; returned: mark if mark > 0, 0 if mark == 0, -1 if mark is negative. */ int mem_stat_mark(mark) int mark; { if (mark < 0) { mem_stat_mark_curr = 0; return -1; /* error */ } else if (mark == 0) { mem_stat_mark_curr = 0; return 0; } mem_stat_mark_curr = mark; mem_stat_mark_many++; return mark; } /* deallocate static variables; Input argument: mark - a positive number denoting the mark; Returned: -1 if mark < 0 (error); 0 if mark == 0; */ int mem_stat_free_list(mark,list) int mark,list; { int i,j; int (*free_fn)(void*); if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS || mem_connect[list].free_funcs == NULL ) return -1; if (mark < 0) { mem_stat_mark_curr = 0; return -1; } else if (mark == 0) { mem_stat_mark_curr = 0; return 0; } if (mem_stat_mark_many <= 0) { warning(WARN_NO_MARK,"mem_stat_free"); return -1; } /* deallocate the marked variables */ for (i=0; i < mem_hash_idx_end; i++) { j = mem_hash_idx[i]; if (j == 0) continue; else { j--; if (mem_stat_var[j].mark == mark) { free_fn = mem_connect[list].free_funcs[mem_stat_var[j].type]; if ( free_fn != NULL ) (*free_fn)(*mem_stat_var[j].var); else warning(WARN_WRONG_TYPE,"mem_stat_free"); *(mem_stat_var[j].var) = NULL; mem_stat_var[j].var = NULL; mem_stat_var[j].mark = 0; mem_hash_idx[i] = 0; } } } while (mem_hash_idx_end > 0 && mem_hash_idx[mem_hash_idx_end-1] == 0) mem_hash_idx_end--; mem_stat_mark_curr = 0; mem_stat_mark_many--; return 0; } /* only for diagnostic purposes */ void mem_stat_dump(fp,list) FILE *fp; int list; { int i,j,k=1; if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS || mem_connect[list].free_funcs == NULL ) return; fprintf(fp," Array mem_stat_var (list no. %d):\n",list); for (i=0; i < mem_hash_idx_end; i++) { j = mem_hash_idx[i]; if (j == 0) continue; else { j--; fprintf(fp," %d. var = 0x%p, type = %s, mark = %d\n", k,mem_stat_var[j].var, mem_stat_var[j].type < mem_connect[list].ntypes && mem_connect[list].free_funcs[mem_stat_var[j].type] != NULL ? mem_connect[list].type_names[(int)mem_stat_var[j].type] : "???", mem_stat_var[j].mark); k++; } } fprintf(fp,"\n"); } /* query function about the current mark */ #ifdef ANSI_C int mem_stat_show_mark(void) #else int mem_stat_show_mark() #endif { return mem_stat_mark_curr; } /* Varying number of arguments */ #ifdef ANSI_C /* To allocate memory to many arguments. The function should be called: mem_stat_vars(list,type,&v1,&v2,&v3,...,VNULL); where int list,type; void **v1, **v2, **v3,...; The last argument should be VNULL ! type is the type of variables v1,v2,v3,... (of course they must be of the same type) */ int mem_stat_reg_vars(int list,int type,...) { va_list ap; int i=0; void **par; va_start(ap, type); while ( (par = va_arg(ap,void **)) ) { /* NULL ends the list*/ mem_stat_reg_list(par,type,list); i++; } va_end(ap); return i; } #elif VARARGS /* old varargs is used */ /* To allocate memory to many arguments. The function should be called: mem_stat_vars(list,type,&v1,&v2,&v3,...,VNULL); where int list,type; void **v1, **v2, **v3,...; The last argument should be VNULL ! type is the type of variables v1,v2,v3,... (of course they must be of the same type) */ int mem_stat_reg_vars(va_alist) va_dcl { va_list ap; int type,list,i=0; void **par; va_start(ap); list = va_arg(ap,int); type = va_arg(ap,int); while (par = va_arg(ap,void **)) { /* NULL ends the list*/ mem_stat_reg_list(par,type,list); i++; } va_end(ap); return i; } #endif /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* File containing routines for determining Hessenberg factorisations. */ /* Hfactor -- compute Hessenberg factorisation in compact form. -- factorisation performed in situ -- for details of the compact form see QRfactor.c and matrix2.doc */ MAT *Hfactor(A, diag, beta) MAT *A; VEC *diag, *beta; { static VEC *tmp1 = VNULL; int k, limit; if ( ! A || ! diag || ! beta ) error(E_NULL,"Hfactor"); if ( diag->dim < A->m - 1 || beta->dim < A->m - 1 ) error(E_SIZES,"Hfactor"); if ( A->m != A->n ) error(E_SQUARE,"Hfactor"); limit = A->m - 1; tmp1 = v_resize(tmp1,A->m); MEM_STAT_REG(tmp1,TYPE_VEC); for ( k = 0; k < limit; k++ ) { get_col(A,k,tmp1); /* printf("the %d'th column = "); v_output(tmp1); */ hhvec(tmp1,k+1,&beta->ve[k],tmp1,&A->me[k+1][k]); /* diag->ve[k] = tmp1->ve[k+1]; */ v_set_val(diag,k,v_entry(tmp1,k+1)); /* printf("H/h vector = "); v_output(tmp1); */ /* printf("from the %d'th entry\n",k+1); */ /* printf("beta = %g\n",beta->ve[k]); */ /* hhtrcols(A,k+1,k+1,tmp1,beta->ve[k]); */ /* hhtrrows(A,0 ,k+1,tmp1,beta->ve[k]); */ hhtrcols(A,k+1,k+1,tmp1,v_entry(beta,k)); hhtrrows(A,0 ,k+1,tmp1,v_entry(beta,k)); /* printf("A = "); m_output(A); */ } return (A); } /* makeHQ -- construct the Hessenberg orthogonalising matrix Q; -- i.e. Hess M = Q.M.Q' */ MAT *makeHQ(H, diag, beta, Qout) MAT *H, *Qout; VEC *diag, *beta; { int i, j, limit; static VEC *tmp1 = VNULL, *tmp2 = VNULL; if ( H==(MAT *)NULL || diag==(VEC *)NULL || beta==(VEC *)NULL ) error(E_NULL,"makeHQ"); limit = H->m - 1; if ( diag->dim < limit || beta->dim < limit ) error(E_SIZES,"makeHQ"); if ( H->m != H->n ) error(E_SQUARE,"makeHQ"); Qout = m_resize(Qout,H->m,H->m); tmp1 = v_resize(tmp1,H->m); tmp2 = v_resize(tmp2,H->m); MEM_STAT_REG(tmp1,TYPE_VEC); MEM_STAT_REG(tmp2,TYPE_VEC); for ( i = 0; i < H->m; i++ ) { /* tmp1 = i'th basis vector */ for ( j = 0; j < H->m; j++ ) /* tmp1->ve[j] = 0.0; */ v_set_val(tmp1,j,0.0); /* tmp1->ve[i] = 1.0; */ v_set_val(tmp1,i,1.0); /* apply H/h transforms in reverse order */ for ( j = limit-1; j >= 0; j-- ) { get_col(H,j,tmp2); /* tmp2->ve[j+1] = diag->ve[j]; */ v_set_val(tmp2,j+1,v_entry(diag,j)); hhtrvec(tmp2,beta->ve[j],j+1,tmp1,tmp1); } /* insert into Qout */ set_col(Qout,(int)i,tmp1); } return (Qout); } /* makeH -- construct actual Hessenberg matrix */ MAT *makeH(H,Hout) MAT *H, *Hout; { int i, j, limit; if ( H==(MAT *)NULL ) error(E_NULL,"makeH"); if ( H->m != H->n ) error(E_SQUARE,"makeH"); Hout = m_resize(Hout,H->m,H->m); Hout = m_copy(H,Hout); limit = H->m; for ( i = 1; i < limit; i++ ) for ( j = 0; j < i-1; j++ ) /* Hout->me[i][j] = 0.0;*/ m_set_val(Hout,i,j,0.0); return (Hout); } /************************************************************************** ** ** Copyright (C) 1993 David E. Stewart & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* This file contains basic routines which are used by the functions in meschach.a etc. These are the routines that should be modified in order to take full advantage of specialised architectures (pipelining, vector processors etc). */ /* __ip__ -- inner product */ double __ip__(dp1,dp2,len) Real *dp1, *dp2; int len; { #ifdef VUNROLL int len4; Real sum1, sum2, sum3; #endif int i; Real sum; sum = 0.0; if (len < 0 ) error(E_BOUNDS,"__ip__"); #ifdef VUNROLL sum1 = sum2 = sum3 = 0.0; len4 = len / 4; len = len % 4; for ( i = 0; i < len4; i++ ) { sum += dp1[4*i]*dp2[4*i]; sum1 += dp1[4*i+1]*dp2[4*i+1]; sum2 += dp1[4*i+2]*dp2[4*i+2]; sum3 += dp1[4*i+3]*dp2[4*i+3]; } sum += sum1 + sum2 + sum3; dp1 += 4*len4; dp2 += 4*len4; #endif for ( i = 0; i < len; i++ ) sum += dp1[i]*dp2[i]; return sum; } /* __mltadd__ -- scalar multiply and add c.f. v_mltadd() */ void __mltadd__(dp1,dp2,s,len) Real *dp1, *dp2; double s; int len; { int i; #ifdef VUNROLL int len4; len4 = len / 4; len = len % 4; for ( i = 0; i < len4; i++ ) { dp1[4*i] += s*dp2[4*i]; dp1[4*i+1] += s*dp2[4*i+1]; dp1[4*i+2] += s*dp2[4*i+2]; dp1[4*i+3] += s*dp2[4*i+3]; } dp1 += 4*len4; dp2 += 4*len4; #endif for ( i = 0; i < len; i++ ) dp1[i] += s*dp2[i]; } /* __smlt__ scalar multiply array c.f. sv_mlt() */ void __smlt__(dp,s,out,len) Real *dp, *out; double s; int len; { int i; for ( i = 0; i < len; i++ ) out[i] = s*dp[i]; } /* __add__ -- add arrays c.f. v_add() */ void __add__(dp1,dp2,out,len) Real *dp1, *dp2, *out; int len; { int i; for ( i = 0; i < len; i++ ) out[i] = dp1[i] + dp2[i]; } /* __sub__ -- subtract arrays c.f. v_sub() */ void __sub__(dp1,dp2,out,len) Real *dp1, *dp2, *out; int len; { int i; for ( i = 0; i < len; i++ ) out[i] = dp1[i] - dp2[i]; } /* __zero__ -- zeros an array of floating point numbers */ void __zero__(dp,len) Real *dp; int len; { #ifdef CHAR0ISDBL0 /* if a floating point zero is equivalent to a string of nulls */ MEM_ZERO((char *)dp,len*sizeof(Real)); #else /* else, need to zero the array entry by entry */ int i; for ( i = 0; i < len; i++ ) dp[i] = 0.0; #endif } /************************************************************************** ** ** Copyright (C) 1993 David E. Stewart & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* File with basic error-handling operations Based on previous version on Zilog System 8000 setret() etc. Ported to Pyramid 9810 late 1987 */ #ifdef SYSV /* AT&T System V */ #include #else /* something else -- assume BSD or ANSI C */ #include #endif #define FALSE 0 #define TRUE 1 #define EF_EXIT 0 #define EF_ABORT 1 #define EF_JUMP 2 #define EF_SILENT 3 /* The only error caught in this file! */ #define E_SIGNAL 16 static char *err_mesg[] = { "unknown error", /* 0 */ "sizes of objects don't match", /* 1 */ "index out of bounds", /* 2 */ "can't allocate memory", /* 3 */ "singular matrix", /* 4 */ "matrix not positive definite", /* 5 */ "incorrect format input", /* 6 */ "bad input file/device", /* 7 */ "NULL objects passed", /* 8 */ "matrix not square", /* 9 */ "object out of range", /* 10 */ "can't do operation in situ for non-square matrix", /* 11 */ "can't do operation in situ", /* 12 */ "excessive number of iterations", /* 13 */ "convergence criterion failed", /* 14 */ "bad starting value", /* 15 */ "floating exception", /* 16 */ "internal inconsistency (data structure)",/* 17 */ "unexpected end-of-file", /* 18 */ "shared vectors (cannot release them)", /* 19 */ "negative argument", /* 20 */ "cannot overwrite object", /* 21 */ "breakdown in iterative method" /* 22 */ }; #define MAXERR (sizeof(err_mesg)/sizeof(char *)) static char *warn_mesg[] = { "unknown warning", /* 0 */ "wrong type number (use macro TYPE_*)", /* 1 */ "no corresponding mem_stat_mark", /* 2 */ "computed norm of a residual is less than 0", /* 3 */ "resizing a shared vector" /* 4 */ }; #define MAXWARN (sizeof(warn_mesg)/sizeof(char *)) #define MAX_ERRS 100 jmp_buf restart; /* array of pointers to lists of errors */ typedef struct { char **listp; /* pointer to a list of errors */ unsigned len; /* length of the list */ unsigned warn; /* =FALSE - errors, =TRUE - warnings */ } Err_list; static Err_list err_list[ERR_LIST_MAX_LEN] = { {err_mesg,MAXERR,FALSE}, /* basic errors list */ {warn_mesg,MAXWARN,TRUE} /* basic warnings list */ }; static int err_list_end = 2; /* number of elements in err_list */ /* attach a new list of errors pointed by err_ptr or change a previous one; list_len is the number of elements in the list; list_num is the list number; warn == FALSE - errors (stop the program), warn == TRUE - warnings (continue the program); Note: lists numbered 0 and 1 are attached automatically, you do not need to do it */ int err_list_attach(list_num, list_len,err_ptr,warn) int list_num, list_len, warn; char **err_ptr; { if (list_num < 0 || list_len <= 0 || err_ptr == (char **)NULL) return -1; if (list_num >= ERR_LIST_MAX_LEN) { fprintf(stderr,"\n file \"%s\": %s %s\n", "err.c","increase the value of ERR_LIST_MAX_LEN", "in matrix.h and zmatdef.h"); if ( ! isatty(fileno(stdout)) ) fprintf(stderr,"\n file \"%s\": %s %s\n", "err.c","increase the value of ERR_LIST_MAX_LEN", "in matrix.h and zmatdef.h"); printf("Exiting program\n"); exit(0); } if (err_list[list_num].listp != (char **)NULL && err_list[list_num].listp != err_ptr) free((char *)err_list[list_num].listp); err_list[list_num].listp = err_ptr; err_list[list_num].len = list_len; err_list[list_num].warn = warn; err_list_end = list_num+1; return list_num; } /* release the error list numbered list_num */ int err_list_free(list_num) int list_num; { if (list_num < 0 || list_num >= err_list_end) return -1; if (err_list[list_num].listp != (char **)NULL) { err_list[list_num].listp = (char **)NULL; err_list[list_num].len = 0; err_list[list_num].warn = 0; } return 0; } /* check if list_num is attached; return FALSE if not; return TRUE if yes */ int err_is_list_attached(list_num) int list_num; { if (list_num < 0 || list_num >= err_list_end) return FALSE; if (err_list[list_num].listp != (char **)NULL) return TRUE; return FALSE; } /* other local variables */ static int err_flag = EF_EXIT, num_errs = 0, cnt_errs = 1; /* set_err_flag -- sets err_flag -- returns old err_flag */ int set_err_flag(flag) int flag; { int tmp; tmp = err_flag; err_flag = flag; return tmp; } /* count_errs -- sets cnt_errs (TRUE/FALSE) & returns old value */ int count_errs(flag) int flag; { int tmp; tmp = cnt_errs; cnt_errs = flag; return tmp; } /* ev_err -- reports error (err_num) in file "file" at line "line_num" and returns to user error handler; list_num is an error list number (0 is the basic list pointed by err_mesg, 1 is the basic list of warnings) */ int ev_err(file,err_num,line_num,fn_name,list_num) char *file, *fn_name; int err_num, line_num,list_num; { int num; if ( err_num < 0 ) err_num = 0; if (list_num < 0 || list_num >= err_list_end || err_list[list_num].listp == (char **)NULL) { fprintf(stderr, "\n Not (properly) attached list of errors: list_num = %d\n", list_num); fprintf(stderr," Call \"err_list_attach\" in your program\n"); if ( ! isatty(fileno(stdout)) ) { fprintf(stderr, "\n Not (properly) attached list of errors: list_num = %d\n", list_num); fprintf(stderr," Call \"err_list_attach\" in your program\n"); } printf("\nExiting program\n"); exit(0); } num = err_num; if ( num >= err_list[list_num].len ) num = 0; if ( cnt_errs && ++num_errs >= MAX_ERRS ) /* too many errors */ { fprintf(stderr,"\n\"%s\", line %d: %s in function %s()\n", file,line_num,err_list[list_num].listp[num], isascii( *fn_name) ? fn_name : "???"); if ( ! isatty(fileno(stdout)) ) fprintf(stdout,"\n\"%s\", line %d: %s in function %s()\n", file,line_num,err_list[list_num].listp[num], isascii(*fn_name) ? fn_name : "???"); printf("Sorry, too many errors: %d\n",num_errs); printf("Exiting program\n"); exit(0); } if ( err_list[list_num].warn ) switch ( err_flag ) { case EF_SILENT: break; default: fprintf(stderr,"\n\"%s\", line %d: %s in function %s()\n\n", file,line_num,err_list[list_num].listp[num], isascii(*fn_name) ? fn_name : "???"); if ( ! isatty(fileno(stdout)) ) fprintf(stdout,"\n\"%s\", line %d: %s in function %s()\n\n", file,line_num,err_list[list_num].listp[num], isascii(*fn_name) ? fn_name : "???"); break; } else switch ( err_flag ) { case EF_SILENT: longjmp(restart,(err_num==0)? -1 : err_num); break; case EF_ABORT: fprintf(stderr,"\n\"%s\", line %d: %s in function %s()\n", file,line_num,err_list[list_num].listp[num], isascii(*fn_name) ? fn_name : "???"); if ( ! isatty(fileno(stdout)) ) fprintf(stdout,"\n\"%s\", line %d: %s in function %s()\n", file,line_num,err_list[list_num].listp[num], isascii(*fn_name) ? fn_name : "???"); abort(); break; case EF_JUMP: fprintf(stderr,"\n\"%s\", line %d: %s in function %s()\n", file,line_num,err_list[list_num].listp[num], isascii(*fn_name) ? fn_name : "???"); if ( ! isatty(fileno(stdout)) ) fprintf(stdout,"\n\"%s\", line %d: %s in function %s()\n", file,line_num,err_list[list_num].listp[num], isascii(*fn_name) ? fn_name : "???"); longjmp(restart,(err_num==0)? -1 : err_num); break; default: fprintf(stderr,"\n\"%s\", line %d: %s in function %s()\n\n", file,line_num,err_list[list_num].listp[num], isascii(*fn_name) ? fn_name : "???"); if ( ! isatty(fileno(stdout)) ) fprintf(stdout,"\n\"%s\", line %d: %s in function %s()\n\n", file,line_num,err_list[list_num].listp[num], isascii(*fn_name) ? fn_name : "???"); break; } /* ensure exit if fall through */ if ( ! err_list[list_num].warn ) exit(0); return 0; } /* float_error -- catches floating arithmetic signals */ static void float_error(int num) { signal(SIGFPE,float_error); /* fprintf(stderr,"SIGFPE: signal #%d\n",num); */ /* fprintf(stderr,"errno = %d\n",errno); */ ev_err("???.c",E_SIGNAL,0,"???",0); } /* catch_signal -- sets up float_error() to catch SIGFPE's */ void catch_FPE(void) { signal(SIGFPE,float_error); } /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* 1.6 matrixio.c 11/25/87 */ /* local variables */ static char line[MAXLINE]; /************************************************************************** Input routines **************************************************************************/ /* skipjunk -- skips white spaces and strings of the form #....\n Here .... is a comment string */ int skipjunk(fp) FILE *fp; { int c; for ( ; ; ) /* forever do... */ { /* skip blanks */ do c = getc(fp); while ( isspace(c) ); /* skip comments (if any) */ if ( c == '#' ) /* yes it is a comment (line) */ while ( (c=getc(fp)) != '\n' ) ; else { ungetc(c,fp); break; } } return 0; } MAT *m_finput(fp,a) FILE *fp; MAT *a; { if ( isatty(fileno(fp)) ) return im_finput(fp,a); else return bm_finput(fp,a); } /* im_finput -- interactive input of matrix */ MAT *im_finput(FILE* fp, MAT* mat) { char c; int i, j, m, n, dynamic; /* dynamic set to TRUE if memory allocated here */ /* get matrix size */ if ( mat != (MAT *)NULL && mat->mnm; n = mat->n; dynamic = FALSE; } else { dynamic = TRUE; do { fprintf(stderr,"Matrix: rows cols:"); if ( fgets(line,MAXLINE,fp)==NULL ) error(E_INPUT,"im_finput"); } while ( sscanf(line,"%u%u",&m,&n)<2 || m>MAXDIM || n>MAXDIM ); mat = m_get(m,n); } /* input elements */ for ( i=0; ime[i][j]); if ( fgets(line,MAXLINE,fp)==NULL ) error(E_INPUT,"im_finput"); if ( (*line == 'b' || *line == 'B') && j > 0 ) { j--; dynamic = FALSE; goto redo2; } if ( (*line == 'f' || *line == 'F') && j < n-1 ) { j++; dynamic = FALSE; goto redo2; } #if REAL == DOUBLE } while ( *line=='\0' || sscanf(line,"%lf",&mat->me[i][j])<1 ); #elif REAL == FLOAT } while ( *line=='\0' || sscanf(line,"%f",&mat->me[i][j])<1 ); #endif fprintf(stderr,"Continue: "); fscanf(fp,"%c",&c); if ( c == 'n' || c == 'N' ) { dynamic = FALSE; goto redo; } if ( (c == 'b' || c == 'B') /* && i > 0 */ ) { if ( i > 0 ) i--; dynamic = FALSE; goto redo; } } return (mat); } /* bm_finput -- batch-file input of matrix */ MAT *bm_finput(fp,mat) FILE *fp; MAT *mat; { int i,j,m,n,dummy; int io_code; /* get dimension */ skipjunk(fp); if ((io_code=fscanf(fp," Matrix: %u by %u",&m,&n)) < 2 || m>MAXDIM || n>MAXDIM ) error(io_code==EOF ? E_EOF : E_FORMAT,"bm_finput"); /* allocate memory if necessary */ if ( mat==(MAT *)NULL ) mat = m_resize(mat,m,n); /* get entries */ for ( i=0; ime[i][j])) < 1 ) #elif REAL == FLOAT if ((io_code=fscanf(fp,"%f",&mat->me[i][j])) < 1 ) #endif error(io_code==EOF ? 7 : 6,"bm_finput"); } return (mat); } PERM *px_finput(fp,px) FILE *fp; PERM *px; { if ( isatty(fileno(fp)) ) return ipx_finput(fp,px); else return bpx_finput(fp,px); } /* ipx_finput -- interactive input of permutation */ PERM *ipx_finput(fp,px) FILE *fp; PERM *px; { int i,j,size,dynamic; /* dynamic set if memory allocated here */ int entry,ok; /* get permutation size */ if ( px!=(PERM *)NULL && px->sizesize; dynamic = FALSE; } else { dynamic = TRUE; do { fprintf(stderr,"Permutation: size: "); if ( fgets(line,MAXLINE,fp)==NULL ) error(E_INPUT,"ipx_finput"); } while ( sscanf(line,"%u",&size)<1 || size>MAXDIM ); px = px_get(size); } /* get entries */ i = 0; while ( i%u new: ", i,px->pe[i]); if ( fgets(line,MAXLINE,fp)==NULL ) error(E_INPUT,"ipx_finput"); if ( (*line == 'b' || *line == 'B') && i > 0 ) { i--; dynamic = FALSE; goto redo; } } while ( *line=='\0' || sscanf(line,"%u",&entry) < 1 ); /* check entry */ ok = (entry < size); for ( j=0; jpe[j]); if ( ok ) { px->pe[i] = entry; i++; } } return (px); } /* bpx_finput -- batch-file input of permutation */ PERM *bpx_finput(fp,px) FILE *fp; PERM *px; { int i,j,size,entry,ok; int io_code; /* get size of permutation */ skipjunk(fp); if ((io_code=fscanf(fp," Permutation: size:%u",&size)) < 1 || size>MAXDIM ) error(io_code==EOF ? 7 : 6,"bpx_finput"); /* allocate memory if necessary */ if ( px==(PERM *)NULL || px->size %u",&entry)) < 1 ) error(io_code==EOF ? 7 : 6,"bpx_finput"); /* check entry */ ok = (entry < size); for ( j=0; jpe[j]); if ( ok ) { px->pe[i] = entry; i++; } else error(E_BOUNDS,"bpx_finput"); } return (px); } VEC *v_finput(fp,x) FILE *fp; VEC *x; { if ( isatty(fileno(fp)) ) return ifin_vec(fp,x); else return bfin_vec(fp,x); } /* ifin_vec -- interactive input of vector */ VEC *ifin_vec(fp,vec) FILE *fp; VEC *vec; { int i,dim,dynamic; /* dynamic set if memory allocated here */ /* get vector dimension */ if ( vec != (VEC *)NULL && vec->dimdim; dynamic = FALSE; } else { dynamic = TRUE; do { fprintf(stderr,"Vector: dim: "); if ( fgets(line,MAXLINE,fp)==NULL ) error(E_INPUT,"ifin_vec"); } while ( sscanf(line,"%u",&dim)<1 || dim>MAXDIM ); vec = v_get(dim); } /* input elements */ for ( i=0; ive[i]); if ( fgets(line,MAXLINE,fp)==NULL ) error(E_INPUT,"ifin_vec"); if ( (*line == 'b' || *line == 'B') && i > 0 ) { i--; dynamic = FALSE; goto redo; } if ( (*line == 'f' || *line == 'F') && i < dim-1 ) { i++; dynamic = FALSE; goto redo; } #if REAL == DOUBLE } while ( *line=='\0' || sscanf(line,"%lf",&vec->ve[i]) < 1 ); #elif REAL == FLOAT } while ( *line=='\0' || sscanf(line,"%f",&vec->ve[i]) < 1 ); #endif return (vec); } /* bfin_vec -- batch-file input of vector */ VEC *bfin_vec(fp,vec) FILE *fp; VEC *vec; { int i,dim; int io_code; /* get dimension */ skipjunk(fp); if ((io_code=fscanf(fp," Vector: dim:%u",&dim)) < 1 || dim>MAXDIM ) error(io_code==EOF ? 7 : 6,"bfin_vec"); /* allocate memory if necessary */ if ( vec==(VEC *)NULL ) vec = v_resize(vec,dim); /* get entries */ skipjunk(fp); for ( i=0; ive[i])) < 1 ) #elif REAL == FLOAT if ((io_code=fscanf(fp,"%f",&vec->ve[i])) < 1 ) #endif error(io_code==EOF ? 7 : 6,"bfin_vec"); return (vec); } /************************************************************************** Output routines **************************************************************************/ static char *format = "%14.9g "; char *setformat(char *f_string) { char *old_f_string; old_f_string = format; if ( f_string != (char *)NULL && *f_string != '\0' ) format = f_string; return old_f_string; } void m_foutput(fp,a) FILE *fp; MAT *a; { int i, j, tmp; if ( a == (MAT *)NULL ) { fprintf(fp,"Matrix: NULL\n"); return; } fprintf(fp,"Matrix: %d by %d\n",a->m,a->n); if ( a->me == (Real **)NULL ) { fprintf(fp,"NULL\n"); return; } for ( i=0; im; i++ ) /* for each row... */ { fprintf(fp,"row %u: ",i); for ( j=0, tmp=2; jn; j++, tmp++ ) { /* for each col in row... */ fprintf(fp,format,a->me[i][j]); if ( ! (tmp % 5) ) putc('\n',fp); } if ( tmp % 5 != 1 ) putc('\n',fp); } } void px_foutput(fp,px) FILE *fp; PERM *px; { int i; if ( px == (PERM *)NULL ) { fprintf(fp,"Permutation: NULL\n"); return; } fprintf(fp,"Permutation: size: %u\n",px->size); if ( px->pe == (int *)NULL ) { fprintf(fp,"NULL\n"); return; } for ( i=0; isize; i++ ) if ( ! (i % 8) && i != 0 ) fprintf(fp,"\n %u->%u ",i,px->pe[i]); else fprintf(fp,"%u->%u ",i,px->pe[i]); fprintf(fp,"\n"); } void v_foutput(fp,x) FILE *fp; VEC *x; { int i, tmp; if ( x == (VEC *)NULL ) { fprintf(fp,"Vector: NULL\n"); return; } fprintf(fp,"Vector: dim: %d\n",x->dim); if ( x->ve == (Real *)NULL ) { fprintf(fp,"NULL\n"); return; } for ( i=0, tmp=0; idim; i++, tmp++ ) { fprintf(fp,format,x->ve[i]); if ( tmp % 5 == 4 ) putc('\n',fp); } if ( tmp % 5 != 0 ) putc('\n',fp); } void m_dump(fp,a) FILE *fp; MAT *a; { int i, j, tmp; if ( a == (MAT *)NULL ) { fprintf(fp,"Matrix: NULL\n"); return; } fprintf(fp,"Matrix: %d by %d @ 0x%lx\n",a->m,a->n,(long)a); fprintf(fp,"\tmax_m = %d, max_n = %d, max_size = %d\n", a->max_m, a->max_n, a->max_size); if ( a->me == (Real **)NULL ) { fprintf(fp,"NULL\n"); return; } fprintf(fp,"a->me @ 0x%lx\n",(long)(a->me)); fprintf(fp,"a->base @ 0x%lx\n",(long)(a->base)); for ( i=0; im; i++ ) /* for each row... */ { fprintf(fp,"row %u: @ 0x%lx ",i,(long)(a->me[i])); for ( j=0, tmp=2; jn; j++, tmp++ ) { /* for each col in row... */ fprintf(fp,format,a->me[i][j]); if ( ! (tmp % 5) ) putc('\n',fp); } if ( tmp % 5 != 1 ) putc('\n',fp); } } void px_dump(fp,px) FILE *fp; PERM *px; { int i; if ( ! px ) { fprintf(fp,"Permutation: NULL\n"); return; } fprintf(fp,"Permutation: size: %u @ 0x%lx\n",px->size,(long)(px)); if ( ! px->pe ) { fprintf(fp,"NULL\n"); return; } fprintf(fp,"px->pe @ 0x%lx\n",(long)(px->pe)); for ( i=0; isize; i++ ) fprintf(fp,"%u->%u ",i,px->pe[i]); fprintf(fp,"\n"); } void v_dump(fp,x) FILE *fp; VEC *x; { int i, tmp; if ( ! x ) { fprintf(fp,"Vector: NULL\n"); return; } fprintf(fp,"Vector: dim: %d @ 0x%lx\n",x->dim,(long)(x)); if ( ! x->ve ) { fprintf(fp,"NULL\n"); return; } fprintf(fp,"x->ve @ 0x%lx\n",(long)(x->ve)); for ( i=0, tmp=0; idim; i++, tmp++ ) { fprintf(fp,format,x->ve[i]); if ( tmp % 5 == 4 ) putc('\n',fp); } if ( tmp % 5 != 0 ) putc('\n',fp); } /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* meminfo.c revised 22/11/93 */ /* contains basic functions, types and arrays to keep track of memory allocation/deallocation */ /* this array is defined further in this file */ extern MEM_CONNECT mem_connect[MEM_CONNECT_MAX_LISTS]; /* attach a new list of types */ int mem_attach_list(list, ntypes, type_names, free_funcs, info_sum) int list,ntypes; /* number of a list and number of types there */ char *type_names[]; /* list of names of types */ int (*free_funcs[])(void*); /* list of releasing functions */ MEM_ARRAY info_sum[]; /* local table */ { if (list < 0 || list >= MEM_CONNECT_MAX_LISTS) return -1; if (type_names == NULL || free_funcs == NULL || info_sum == NULL || ntypes < 0) return -1; /* if a list exists do not overwrite */ if ( mem_connect[list].ntypes != 0 ) error(E_OVERWRITE,"mem_attach_list"); mem_connect[list].ntypes = ntypes; mem_connect[list].type_names = type_names; mem_connect[list].free_funcs = free_funcs; mem_connect[list].info_sum = info_sum; return 0; } /* release a list of types */ int mem_free_vars(int list) { if (list < 0 || list >= MEM_CONNECT_MAX_LISTS) return -1; mem_connect[list].ntypes = 0; mem_connect[list].type_names = NULL; mem_connect[list].free_funcs = NULL; mem_connect[list].info_sum = NULL; return 0; } /* check if list is attached */ int mem_is_list_attached(int list) { if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS ) return FALSE; if ( mem_connect[list].type_names != NULL && mem_connect[list].free_funcs != NULL && mem_connect[list].info_sum != NULL) return TRUE; else return FALSE; } /* to print out the contents of mem_connect[list] */ void mem_dump_list(FILE* fp, int list) { int i; MEM_CONNECT *mlist; if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS ) return; mlist = &mem_connect[list]; fprintf(fp," %15s[%d]:\n","CONTENTS OF mem_connect",list); fprintf(fp," %-7s %-12s %-9s %s\n", "name of", "alloc.", "# alloc.", "address" ); fprintf(fp," %-7s %-12s %-9s %s\n", " type", "bytes", "variables", "of *_free()" ); for (i=0; i < mlist->ntypes; i++) fprintf(fp," %-7s %-12ld %-9d %p\n", mlist->type_names[i], mlist->info_sum[i].bytes, mlist->info_sum[i].numvar, mlist->free_funcs[i] ); fprintf(fp,"\n"); } /*=============================================================*/ /* local variables */ static int mem_switched_on = MEM_SWITCH_ON_DEF; /* on/off */ /* switch on/off memory info */ int mem_info_on(int sw) { int old = mem_switched_on; mem_switched_on = sw; return old; } #ifdef ANSI_C int mem_info_is_on(void) #else int mem_info_is_on() #endif { return mem_switched_on; } /* information about allocated memory */ /* return the number of allocated bytes for type 'type' */ long mem_info_bytes(int type, int list) { if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS ) return 0l; if ( !mem_switched_on || type < 0 || type >= mem_connect[list].ntypes || mem_connect[list].free_funcs[type] == NULL ) return 0l; return mem_connect[list].info_sum[type].bytes; } /* return the number of allocated variables for type 'type' */ int mem_info_numvar(int type, int list) { if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS ) return 0l; if ( !mem_switched_on || type < 0 || type >= mem_connect[list].ntypes || mem_connect[list].free_funcs[type] == NULL ) return 0l; return mem_connect[list].info_sum[type].numvar; } /* print out memory info to the file fp */ void mem_info_file(FILE* fp, int list) { unsigned int type; long t = 0l, d; int n = 0, nt = 0; MEM_CONNECT *mlist; if (!mem_switched_on) return; if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS ) return; if (list == 0) fprintf(fp," MEMORY INFORMATION (standard types):\n"); else fprintf(fp," MEMORY INFORMATION (list no. %d):\n",list); mlist = &mem_connect[list]; for (type=0; type < mlist->ntypes; type++) { if (mlist->type_names[type] == NULL ) continue; d = mlist->info_sum[type].bytes; t += d; n = mlist->info_sum[type].numvar; nt += n; fprintf(fp," type %-7s %10ld alloc. byte%c %6d alloc. variable%c\n", mlist->type_names[type], d, (d!=1 ? 's' : ' '), n, (n!=1 ? 's' : ' ')); } fprintf(fp," %-12s %10ld alloc. byte%c %6d alloc. variable%c\n\n", "total:",t, (t!=1 ? 's' : ' '), nt, (nt!=1 ? 's' : ' ')); } /* function for memory information */ /* mem_bytes_list Arguments: type - the number of type; old_size - old size of allocated memory (in bytes); new_size - new size of allocated memory (in bytes); list - list of types */ void mem_bytes_list(int type, int old_size, int new_size, int list) { MEM_CONNECT *mlist; if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS ) return; mlist = &mem_connect[list]; if ( type < 0 || type >= mlist->ntypes || mlist->free_funcs[type] == NULL ) return; if ( old_size < 0 || new_size < 0 ) error(E_NEG,"mem_bytes_list"); mlist->info_sum[type].bytes += new_size - old_size; /* check if the number of bytes is non-negative */ if ( old_size > 0 ) { if (mlist->info_sum[type].bytes < 0) { fprintf(stderr, "\n WARNING !! memory info: allocated memory is less than 0\n"); fprintf(stderr,"\t TYPE %s \n\n", mlist->type_names[type]); if ( !isatty(fileno(stdout)) ) { fprintf(stdout, "\n WARNING !! memory info: allocated memory is less than 0\n"); fprintf(stdout,"\t TYPE %s \n\n", mlist->type_names[type]); } } } } /* mem_numvar_list Arguments: type - the number of type; num - # of variables allocated (> 0) or deallocated ( < 0) list - list of types */ void mem_numvar_list(int type, int num, int list) { MEM_CONNECT *mlist; if ( list < 0 || list >= MEM_CONNECT_MAX_LISTS ) return; mlist = &mem_connect[list]; if ( type < 0 || type >= mlist->ntypes || mlist->free_funcs[type] == NULL ) return; mlist->info_sum[type].numvar += num; /* check if the number of variables is non-negative */ if ( num < 0 ) { if (mlist->info_sum[type].numvar < 0) { fprintf(stderr, "\n WARNING !! memory info: allocated # of variables is less than 0\n"); fprintf(stderr,"\t TYPE %s \n\n", mlist->type_names[type]); if ( !isatty(fileno(stdout)) ) { fprintf(stdout, "\n WARNING !! memory info: allocated # of variables is less than 0\n"); fprintf(stdout,"\t TYPE %s \n\n", mlist->type_names[type]); } } } } /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* This is a file of routines for zero-ing, and initialising vectors, matrices and permutations. This is to be included in the matrix.a library */ /* v_zero -- zero the vector x */ VEC *v_zero(VEC *x) { if ( x == VNULL ) error(E_NULL,"v_zero"); __zero__(x->ve,x->dim); /* for ( i = 0; i < x->dim; i++ ) x->ve[i] = 0.0; */ return x; } /* iv_zero -- zero the vector ix */ IVEC *iv_zero(ix) IVEC *ix; { int i; if ( ix == IVNULL ) error(E_NULL,"iv_zero"); for ( i = 0; i < ix->dim; i++ ) ix->ive[i] = 0; return ix; } /* m_zero -- zero the matrix A */ MAT *m_zero(A) MAT *A; { int i, A_m, A_n; Real **A_me; if ( A == MNULL ) error(E_NULL,"m_zero"); A_m = A->m; A_n = A->n; A_me = A->me; for ( i = 0; i < A_m; i++ ) __zero__(A_me[i],A_n); /* for ( j = 0; j < A_n; j++ ) A_me[i][j] = 0.0; */ return A; } /* mat_id -- set A to being closest to identity matrix as possible -- i.e. A[i][j] == 1 if i == j and 0 otherwise */ MAT *m_ident(A) MAT *A; { int i, size; if ( A == MNULL ) error(E_NULL,"m_ident"); m_zero(A); size = min(A->m,A->n); for ( i = 0; i < size; i++ ) A->me[i][i] = 1.0; return A; } /* px_ident -- set px to identity permutation */ PERM *px_ident(px) PERM *px; { int i, px_size; int *px_pe; if ( px == PNULL ) error(E_NULL,"px_ident"); px_size = px->size; px_pe = px->pe; for ( i = 0; i < px_size; i++ ) px_pe[i] = i; return px; } /* Pseudo random number generator data structures */ /* Knuth's lagged Fibonacci-based generator: See "Seminumerical Algorithms: The Art of Computer Programming" sections 3.2-3.3 */ #ifdef ANSI_C #ifndef LONG_MAX #include #endif #endif #ifdef LONG_MAX #define MODULUS LONG_MAX #else #define MODULUS 1000000000L /* assuming long's at least 32 bits long */ #endif #define MZ 0L static long mrand_list[56]; static int started = FALSE; static int inext = 0, inextp = 31; /* mrand -- pseudo-random number generator */ #ifdef ANSI_C double mrand(void) #else double mrand() #endif { long lval; static Real factor = 1.0/((Real)MODULUS); if ( ! started ) smrand(3127); inext = (inext >= 54) ? 0 : inext+1; inextp = (inextp >= 54) ? 0 : inextp+1; lval = mrand_list[inext]-mrand_list[inextp]; if ( lval < 0L ) lval += MODULUS; mrand_list[inext] = lval; return (double)lval*factor; } /* mrandlist -- fills the array a[] with len random numbers */ void mrandlist(a, len) Real a[]; int len; { int i; long lval; static Real factor = 1.0/((Real)MODULUS); if ( ! started ) smrand(3127); for ( i = 0; i < len; i++ ) { inext = (inext >= 54) ? 0 : inext+1; inextp = (inextp >= 54) ? 0 : inextp+1; lval = mrand_list[inext]-mrand_list[inextp]; if ( lval < 0L ) lval += MODULUS; mrand_list[inext] = lval; a[i] = (Real)lval*factor; } } /* smrand -- set seed for mrand() */ void smrand(seed) int seed; { int i; mrand_list[0] = (123413*seed) % MODULUS; for ( i = 1; i < 55; i++ ) mrand_list[i] = (123413*mrand_list[i-1]) % MODULUS; started = TRUE; /* run mrand() through the list sufficient times to thoroughly randomise the array */ for ( i = 0; i < 55*55; i++ ) mrand(); } #undef MODULUS #undef MZ #undef FAC /* v_rand -- initialises x to be a random vector, components independently & uniformly ditributed between 0 and 1 */ VEC *v_rand(x) VEC *x; { /* int i; */ if ( ! x ) error(E_NULL,"v_rand"); /* for ( i = 0; i < x->dim; i++ ) */ /* x->ve[i] = rand()/((Real)MAX_RAND); */ /* x->ve[i] = mrand(); */ mrandlist(x->ve,x->dim); return x; } /* m_rand -- initialises A to be a random vector, components independently & uniformly distributed between 0 and 1 */ MAT *m_rand(A) MAT *A; { int i /* , j */; if ( ! A ) error(E_NULL,"m_rand"); for ( i = 0; i < A->m; i++ ) /* for ( j = 0; j < A->n; j++ ) */ /* A->me[i][j] = rand()/((Real)MAX_RAND); */ /* A->me[i][j] = mrand(); */ mrandlist(A->me[i],A->n); return A; } /* v_ones -- fills x with one's */ VEC *v_ones(x) VEC *x; { int i; if ( ! x ) error(E_NULL,"v_ones"); for ( i = 0; i < x->dim; i++ ) x->ve[i] = 1.0; return x; } /* m_ones -- fills matrix with one's */ MAT *m_ones(A) MAT *A; { int i, j; if ( ! A ) error(E_NULL,"m_ones"); for ( i = 0; i < A->m; i++ ) for ( j = 0; j < A->n; j++ ) A->me[i][j] = 1.0; return A; } /* v_count -- initialises x so that x->ve[i] == i */ VEC *v_count(VEC *x) { int i; if ( ! x ) error(E_NULL,"v_count"); for ( i = 0; i < x->dim; i++ ) x->ve[i] = (Real)i; return x; } /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* 1.2 submat.c 11/25/87 */ /* get_col -- gets a specified column of a matrix and retruns it as a vector */ VEC *get_col(mat,col,vec) int col; MAT *mat; VEC *vec; { int i; if ( mat==(MAT *)NULL ) error(E_NULL,"get_col"); if ( col >= mat->n ) error(E_RANGE,"get_col"); if ( col < 0 ) error(E_RANGE,"get_col"); if ( vec==(VEC *)NULL || vec->dimm ) vec = v_resize(vec,mat->m); for ( i=0; im; i++ ) vec->ve[i] = mat->me[i][col]; return (vec); } /* get_row -- gets a specified row of a matrix and retruns it as a vector */ VEC *get_row(mat,row,vec) int row; MAT *mat; VEC *vec; { int i; if ( mat==(MAT *)NULL ) error(E_NULL,"get_row"); if ( row >= mat->m ) error(E_RANGE,"get_row"); if ( row < 0 ) error(E_RANGE,"get_row"); if ( vec==(VEC *)NULL || vec->dimn ) vec = v_resize(vec,mat->n); for ( i=0; in; i++ ) vec->ve[i] = mat->me[row][i]; return (vec); } /* _set_col -- sets column of matrix to values given in vec (in situ) */ MAT *_set_col(mat,col,vec,i0) MAT *mat; VEC *vec; int col,i0; { int i,lim; if ( mat==(MAT *)NULL || vec==(VEC *)NULL ) error(E_NULL,"_set_col"); if ( col >= mat->n ) error(E_RANGE,"_set_col"); lim = min(mat->m,vec->dim); for ( i=i0; ime[i][col] = vec->ve[i]; return (mat); } /* _set_row -- sets row of matrix to values given in vec (in situ) */ MAT *_set_row(mat,row,vec,j0) MAT *mat; VEC *vec; int row,j0; { int j,lim; if ( mat==(MAT *)NULL || vec==(VEC *)NULL ) error(E_NULL,"_set_row"); if ( row >= mat->m ) error(E_RANGE,"_set_row"); lim = min(mat->n,vec->dim); for ( j=j0; jme[row][j] = vec->ve[j]; return (mat); } /* sub_mat -- returns sub-matrix of old which is formed by the rectangle from (row1,col1) to (row2,col2) -- Note: storage is shared so that altering the "new" matrix will alter the "old" matrix */ MAT *sub_mat(old,row1,col1,row2,col2,new) MAT *old,*new; int row1,col1,row2,col2; { int i; if ( old==(MAT *)NULL ) error(E_NULL,"sub_mat"); if ( row1 > row2 || col1 > col2 || row2 >= old->m || col2 >= old->n ) error(E_RANGE,"sub_mat"); if ( new==(MAT *)NULL || new->m < row2-row1+1 ) { new = NEW(MAT); new->me = NEW_A(row2-row1+1,Real *); if ( new==(MAT *)NULL || new->me==(Real **)NULL ) error(E_MEM,"sub_mat"); else if (mem_info_is_on()) { mem_bytes(TYPE_MAT,0,(signed)((sizeof(MAT)+ (row2-row1+1)*sizeof(Real *)))); } } new->m = row2-row1+1; new->n = col2-col1+1; new->base = (Real *)NULL; for ( i=0; i < new->m; i++ ) new->me[i] = (old->me[i+row1]) + col1; return (new); } /* sub_vec -- returns sub-vector which is formed by the elements i1 to i2 -- as for sub_mat, storage is shared */ VEC *sub_vec(old,i1,i2,new) VEC *old, *new; int i1, i2; { if ( old == (VEC *)NULL ) error(E_NULL,"sub_vec"); if ( i2 < i1 || i1 < 0 || old->dim < i2 ) error(E_RANGE,"sub_vec"); if ( new == (VEC *)NULL ) new = NEW(VEC); if ( new == (VEC *)NULL ) error(E_MEM,"sub_vec"); else if (mem_info_is_on()) { mem_bytes(TYPE_VEC,0,sizeof(VEC)); } new->dim = i2 - i1 + 1; new->ve = &(old->ve[i1]); return new; } /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* Files for matrix computations Householder transformation file. Contains routines for calculating householder transformations, applying them to vectors and matrices by both row & column. */ /* hsehldr.c 1.3 10/8/87 */ /* hhvec -- calulates Householder vector to eliminate all entries after the i0 entry of the vector vec. It is returned as out. May be in-situ */ VEC *hhvec(vec,i0,beta,out,newval) VEC *vec,*out; int i0; Real *beta,*newval; { Real norm; if ( i0 < 0 || out->dim < i0 ) error(E_RANGE,"hhvec"); out = _v_copy(vec,out,i0); norm = sqrt(_in_prod(out,out,i0)); if ( norm <= 0.0 ) { *beta = 0.0; return (out); } *beta = 1.0/(norm * (norm+fabs(out->ve[i0]))); if ( out->ve[i0] > 0.0 ) *newval = -norm; else *newval = norm; out->ve[i0] -= *newval; return (out); } /* hhtrvec -- apply Householder transformation to vector -- may be in-situ */ VEC *hhtrvec(hh,beta,i0,in,out) VEC *hh,*in,*out; /* hh = Householder vector */ int i0; double beta; { Real scale; /* int i; */ if ( hh==(VEC *)NULL || in==(VEC *)NULL ) error(E_NULL,"hhtrvec"); if ( in->dim != hh->dim ) error(E_SIZES,"hhtrvec"); if ( i0 < 0 || in->dim < i0 ) error(E_BOUNDS,"hhtrvec"); scale = beta*_in_prod(hh,in,i0); out = v_copy(in,out); __mltadd__(&(out->ve[i0]),&(hh->ve[i0]),-scale,(in->dim-i0)); /************************************************************ for ( i=i0; idim; i++ ) out->ve[i] = in->ve[i] - scale*hh->ve[i]; ************************************************************/ return (out); } /* hhtrrows -- transform a matrix by a Householder vector by rows starting at row i0 from column j0 -- in-situ */ MAT *hhtrrows(M,i0,j0,hh,beta) MAT *M; int i0, j0; VEC *hh; double beta; { Real ip, scale; int i /*, j */; if ( M==(MAT *)NULL || hh==(VEC *)NULL ) error(E_NULL,"hhtrrows"); if ( M->n != hh->dim ) error(E_RANGE,"hhtrrows"); if ( i0 < 0 || M->m < i0 || j0 < 0 || M->n < j0 ) error(E_BOUNDS,"hhtrrows"); if ( beta == 0.0 ) return (M); /* for each row ... */ for ( i = i0; i < M->m; i++ ) { /* compute inner product */ ip = __ip__(&(M->me[i][j0]),&(hh->ve[j0]),(int)(M->n-j0)); /************************************************** ip = 0.0; for ( j = j0; j < M->n; j++ ) ip += M->me[i][j]*hh->ve[j]; **************************************************/ scale = beta*ip; if ( scale == 0.0 ) continue; /* do operation */ __mltadd__(&(M->me[i][j0]),&(hh->ve[j0]),-scale, (M->n-j0)); /************************************************** for ( j = j0; j < M->n; j++ ) M->me[i][j] -= scale*hh->ve[j]; **************************************************/ } return (M); } /* hhtrcols -- transform a matrix by a Householder vector by columns starting at row i0 from column j0 -- in-situ */ MAT *hhtrcols(M,i0,j0,hh,beta) MAT *M; int i0, j0; VEC *hh; double beta; { /* Real ip, scale; */ int i /*, k */; static VEC *w = VNULL; if ( M==(MAT *)NULL || hh==(VEC *)NULL ) error(E_NULL,"hhtrcols"); if ( M->m != hh->dim ) error(E_SIZES,"hhtrcols"); if ( i0 < 0 || M->m < i0 || j0 < 0 || M->n < j0 ) error(E_BOUNDS,"hhtrcols"); if ( beta == 0.0 ) return (M); w = v_resize(w,M->n); MEM_STAT_REG(w,TYPE_VEC); v_zero(w); for ( i = i0; i < M->m; i++ ) if ( hh->ve[i] != 0.0 ) __mltadd__(&(w->ve[j0]),&(M->me[i][j0]),hh->ve[i], (M->n-j0)); for ( i = i0; i < M->m; i++ ) if ( hh->ve[i] != 0.0 ) __mltadd__(&(M->me[i][j0]),&(w->ve[j0]),-beta*hh->ve[i], (M->n-j0)); return (M); } /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* vecop.c 1.3 8/18/87 */ /* _in_prod -- inner product of two vectors from i0 downwards */ double _in_prod(a,b,i0) VEC *a,*b; int i0; { int limit; /* Real *a_v, *b_v; */ /* register Real sum; */ if ( a==(VEC *)NULL || b==(VEC *)NULL ) error(E_NULL,"_in_prod"); limit = min(a->dim,b->dim); if ( i0 < 0 || limit < i0 ) error(E_BOUNDS,"_in_prod"); return __ip__(&(a->ve[i0]),&(b->ve[i0]),(int)(limit-i0)); /***************************************** a_v = &(a->ve[i0]); b_v = &(b->ve[i0]); for ( i=i0; idim != vector->dim ) out = v_resize(out,vector->dim); if ( scalar == 0.0 ) return v_zero(out); if ( scalar == 1.0 ) return v_copy(vector,out); __smlt__(vector->ve,(double)scalar,out->ve,(vector->dim)); /************************************************** dim = vector->dim; out_ve = out->ve; vec_ve = vector->ve; for ( i=0; ive[i] = scalar*vector->ve[i]; (*out_ve++) = scalar*(*vec_ve++); **************************************************/ return (out); } /* v_add -- vector addition -- may be in-situ */ VEC *v_add(vec1,vec2,out) VEC *vec1,*vec2,*out; { int dim; /* Real *out_ve, *vec1_ve, *vec2_ve; */ if ( vec1==(VEC *)NULL || vec2==(VEC *)NULL ) error(E_NULL,"v_add"); if ( vec1->dim != vec2->dim ) error(E_SIZES,"v_add"); if ( out==(VEC *)NULL || out->dim != vec1->dim ) out = v_resize(out,vec1->dim); dim = vec1->dim; __add__(vec1->ve,vec2->ve,out->ve,dim); /************************************************************ out_ve = out->ve; vec1_ve = vec1->ve; vec2_ve = vec2->ve; for ( i=0; ive[i] = vec1->ve[i]+vec2->ve[i]; (*out_ve++) = (*vec1_ve++) + (*vec2_ve++); ************************************************************/ return (out); } /* v_mltadd -- scalar/vector multiplication and addition -- out = v1 + scale.v2 */ VEC *v_mltadd(v1,v2,scale,out) VEC *v1,*v2,*out; double scale; { /* register int dim, i; */ /* Real *out_ve, *v1_ve, *v2_ve; */ if ( v1==(VEC *)NULL || v2==(VEC *)NULL ) error(E_NULL,"v_mltadd"); if ( v1->dim != v2->dim ) error(E_SIZES,"v_mltadd"); if ( scale == 0.0 ) return v_copy(v1,out); if ( scale == 1.0 ) return v_add(v1,v2,out); if ( v2 != out ) { tracecatch(out = v_copy(v1,out),"v_mltadd"); /* dim = v1->dim; */ __mltadd__(out->ve,v2->ve,scale,(v1->dim)); } else { tracecatch(out = sv_mlt(scale,v2,out),"v_mltadd"); out = v_add(v1,out,out); } /************************************************************ out_ve = out->ve; v1_ve = v1->ve; v2_ve = v2->ve; for ( i=0; i < dim ; i++ ) out->ve[i] = v1->ve[i] + scale*v2->ve[i]; (*out_ve++) = (*v1_ve++) + scale*(*v2_ve++); ************************************************************/ return (out); } /* v_sub -- vector subtraction -- may be in-situ */ VEC *v_sub(vec1,vec2,out) VEC *vec1,*vec2,*out; { /* int i, dim; */ /* Real *out_ve, *vec1_ve, *vec2_ve; */ if ( vec1==(VEC *)NULL || vec2==(VEC *)NULL ) error(E_NULL,"v_sub"); if ( vec1->dim != vec2->dim ) error(E_SIZES,"v_sub"); if ( out==(VEC *)NULL || out->dim != vec1->dim ) out = v_resize(out,vec1->dim); __sub__(vec1->ve,vec2->ve,out->ve,(vec1->dim)); /************************************************************ dim = vec1->dim; out_ve = out->ve; vec1_ve = vec1->ve; vec2_ve = vec2->ve; for ( i=0; ive[i] = vec1->ve[i]-vec2->ve[i]; (*out_ve++) = (*vec1_ve++) - (*vec2_ve++); ************************************************************/ return (out); } /* v_map -- maps function f over components of x: out[i] = f(x[i]) -- _v_map sets out[i] = f(params,x[i]) */ VEC *v_map(f,x,out) #ifdef PROTOTYPES_IN_STRUCT double (*f)(double); #else double (*f)(); #endif VEC *x, *out; { Real *x_ve, *out_ve; int i, dim; if ( ! x || ! f ) error(E_NULL,"v_map"); if ( ! out || out->dim != x->dim ) out = v_resize(out,x->dim); dim = x->dim; x_ve = x->ve; out_ve = out->ve; for ( i = 0; i < dim; i++ ) *out_ve++ = (*f)(*x_ve++); return out; } VEC *_v_map(f,params,x,out) #ifdef PROTOTYPES_IN_STRUCT double (*f)(void *,double); #else double (*f)(); #endif VEC *x, *out; void *params; { Real *x_ve, *out_ve; int i, dim; if ( ! x || ! f ) error(E_NULL,"_v_map"); if ( ! out || out->dim != x->dim ) out = v_resize(out,x->dim); dim = x->dim; x_ve = x->ve; out_ve = out->ve; for ( i = 0; i < dim; i++ ) *out_ve++ = (*f)(params,*x_ve++); return out; } /* v_lincomb -- returns sum_i a[i].v[i], a[i] real, v[i] vectors */ VEC *v_lincomb(n,v,a,out) int n; /* number of a's and v's */ Real a[]; VEC *v[], *out; { int i; if ( ! a || ! v ) error(E_NULL,"v_lincomb"); if ( n <= 0 ) return VNULL; for ( i = 1; i < n; i++ ) if ( out == v[i] ) error(E_INSITU,"v_lincomb"); out = sv_mlt(a[0],v[0],out); for ( i = 1; i < n; i++ ) { if ( ! v[i] ) error(E_NULL,"v_lincomb"); if ( v[i]->dim != out->dim ) error(E_SIZES,"v_lincomb"); out = v_mltadd(out,v[i],a[i],out); } return out; } #ifdef ANSI_C /* v_linlist -- linear combinations taken from a list of arguments; calling: v_linlist(out,v1,a1,v2,a2,...,vn,an,NULL); where vi are vectors (VEC *) and ai are numbers (double) */ VEC *v_linlist(VEC *out,VEC *v1,double a1,...) { va_list ap; VEC *par; double a_par; if ( ! v1 ) return VNULL; va_start(ap, a1); out = sv_mlt(a1,v1,out); while ( (par = va_arg(ap,VEC *)) ) { /* NULL ends the list*/ a_par = va_arg(ap,double); if (a_par == 0.0) continue; if ( out == par ) error(E_INSITU,"v_linlist"); if ( out->dim != par->dim ) error(E_SIZES,"v_linlist"); if (a_par == 1.0) out = v_add(out,par,out); else if (a_par == -1.0) out = v_sub(out,par,out); else out = v_mltadd(out,par,a_par,out); } va_end(ap); return out; } #elif VARARGS /* v_linlist -- linear combinations taken from a list of arguments; calling: v_linlist(out,v1,a1,v2,a2,...,vn,an,NULL); where vi are vectors (VEC *) and ai are numbers (double) */ VEC *v_linlist(va_alist) va_dcl { va_list ap; VEC *par, *out; double a_par; va_start(ap); out = va_arg(ap,VEC *); par = va_arg(ap,VEC *); if ( ! par ) { va_end(ap); return VNULL; } a_par = va_arg(ap,double); out = sv_mlt(a_par,par,out); while (par = va_arg(ap,VEC *)) { /* NULL ends the list*/ a_par = va_arg(ap,double); if (a_par == 0.0) continue; if ( out == par ) error(E_INSITU,"v_linlist"); if ( out->dim != par->dim ) error(E_SIZES,"v_linlist"); if (a_par == 1.0) out = v_add(out,par,out); else if (a_par == -1.0) out = v_sub(out,par,out); else out = v_mltadd(out,par,a_par,out); } va_end(ap); return out; } #endif /* v_star -- computes componentwise (Hadamard) product of x1 and x2 -- result out is returned */ VEC *v_star(x1, x2, out) VEC *x1, *x2, *out; { int i; if ( ! x1 || ! x2 ) error(E_NULL,"v_star"); if ( x1->dim != x2->dim ) error(E_SIZES,"v_star"); out = v_resize(out,x1->dim); for ( i = 0; i < x1->dim; i++ ) out->ve[i] = x1->ve[i] * x2->ve[i]; return out; } /* v_slash -- computes componentwise ratio of x2 and x1 -- out[i] = x2[i] / x1[i] -- if x1[i] == 0 for some i, then raise E_SING error -- result out is returned */ VEC *v_slash(x1, x2, out) VEC *x1, *x2, *out; { int i; Real tmp; if ( ! x1 || ! x2 ) error(E_NULL,"v_slash"); if ( x1->dim != x2->dim ) error(E_SIZES,"v_slash"); out = v_resize(out,x1->dim); for ( i = 0; i < x1->dim; i++ ) { tmp = x1->ve[i]; if ( tmp == 0.0 ) error(E_SING,"v_slash"); out->ve[i] = x2->ve[i] / tmp; } return out; } /* v_min -- computes minimum component of x, which is returned -- also sets min_idx to the index of this minimum */ double v_min(x, min_idx) VEC *x; int *min_idx; { int i, i_min; Real min_val, tmp; if ( ! x ) error(E_NULL,"v_min"); if ( x->dim <= 0 ) error(E_SIZES,"v_min"); i_min = 0; min_val = x->ve[0]; for ( i = 1; i < x->dim; i++ ) { tmp = x->ve[i]; if ( tmp < min_val ) { min_val = tmp; i_min = i; } } if ( min_idx != NULL ) *min_idx = i_min; return min_val; } /* v_max -- computes maximum component of x, which is returned -- also sets max_idx to the index of this maximum */ double v_max(x, max_idx) VEC *x; int *max_idx; { int i, i_max; Real max_val, tmp; if ( ! x ) error(E_NULL,"v_max"); if ( x->dim <= 0 ) error(E_SIZES,"v_max"); i_max = 0; max_val = x->ve[0]; for ( i = 1; i < x->dim; i++ ) { tmp = x->ve[i]; if ( tmp > max_val ) { max_val = tmp; i_max = i; } } if ( max_idx != NULL ) *max_idx = i_max; return max_val; } #define MAX_STACK 60 /* v_sort -- sorts vector x, and generates permutation that gives the order of the components; x = [1.3, 3.7, 0.5] -> [0.5, 1.3, 3.7] and the permutation is order = [2, 0, 1]. -- if order is NULL on entry then it is ignored -- the sorted vector x is returned */ VEC *v_sort(x, order) VEC *x; PERM *order; { Real *x_ve, tmp, v; /* int *order_pe; */ int dim, i, j, l, r, tmp_i; int stack[MAX_STACK], sp; if ( ! x ) error(E_NULL,"v_sort"); if ( order != PNULL && order->size != x->dim ) order = px_resize(order, x->dim); x_ve = x->ve; dim = x->dim; if ( order != PNULL ) px_ident(order); if ( dim <= 1 ) return x; /* using quicksort algorithm in Sedgewick, "Algorithms in C", Ch. 9, pp. 118--122 (1990) */ sp = 0; l = 0; r = dim-1; v = x_ve[0]; for ( ; ; ) { while ( r > l ) { /* "i = partition(x_ve,l,r);" */ v = x_ve[r]; i = l-1; j = r; for ( ; ; ) { while ( x_ve[++i] < v ) ; while ( x_ve[--j] > v ) ; if ( i >= j ) break; tmp = x_ve[i]; x_ve[i] = x_ve[j]; x_ve[j] = tmp; if ( order != PNULL ) { tmp_i = order->pe[i]; order->pe[i] = order->pe[j]; order->pe[j] = tmp_i; } } tmp = x_ve[i]; x_ve[i] = x_ve[r]; x_ve[r] = tmp; if ( order != PNULL ) { tmp_i = order->pe[i]; order->pe[i] = order->pe[r]; order->pe[r] = tmp_i; } if ( i-l > r-i ) { stack[sp++] = l; stack[sp++] = i-1; l = i+1; } else { stack[sp++] = i+1; stack[sp++] = r; r = i-1; } } /* recursion elimination */ if ( sp == 0 ) break; r = stack[--sp]; l = stack[--sp]; } return x; } /* v_sum -- returns sum of entries of a vector */ double v_sum(x) VEC *x; { int i; Real sum; if ( ! x ) error(E_NULL,"v_sum"); sum = 0.0; for ( i = 0; i < x->dim; i++ ) sum += x->ve[i]; return sum; } /* v_conv -- computes convolution product of two vectors */ VEC *v_conv(VEC *x1, VEC *x2, VEC *out) { int i; if ( ! x1 || ! x2 ) error(E_NULL,"v_conv"); if ( x1 == out || x2 == out ) error(E_INSITU,"v_conv"); if ( x1->dim == 0 || x2->dim == 0 ) return out = v_resize(out,0); out = v_resize(out,x1->dim + x2->dim - 1); v_zero(out); for ( i = 0; i < x1->dim; i++ ) __mltadd__(&(out->ve[i]),x2->ve,x1->ve[i],x2->dim); return out; } /* v_pconv -- computes a periodic convolution product -- the period is the dimension of x2 */ VEC *v_pconv(VEC *x1, VEC *x2, VEC *out) { int i; if ( ! x1 || ! x2 ) error(E_NULL,"v_pconv"); if ( x1 == out || x2 == out ) error(E_INSITU,"v_pconv"); out = v_resize(out,x2->dim); if ( x2->dim == 0 ) return out; v_zero(out); for ( i = 0; i < x1->dim; i++ ) { __mltadd__(&(out->ve[i]),x2->ve,x1->ve[i],x2->dim - i); if ( i > 0 ) __mltadd__(out->ve,&(x2->ve[x2->dim - i]),x1->ve[i],i); } return out; } /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* Band matrix factorisation routines */ /* bdfactor.c 18/11/93 */ /* generate band matrix for a matrix with n columns, lb subdiagonals and ub superdiagonals; Way of saving a band of a matrix: first we save subdiagonals (from 0 to lb-1); then main diagonal (in the lb row) and then superdiagonals (from lb+1 to lb+ub) in such a way that the elements which were previously in one column are now also in one column */ BAND *bd_get(lb,ub,n) int lb, ub, n; { BAND *A; if (lb < 0 || ub < 0 || n <= 0) error(E_NEG,"bd_get"); if ((A = NEW(BAND)) == (BAND *)NULL) error(E_MEM,"bd_get"); else if (mem_info_is_on()) { mem_bytes(TYPE_BAND,0,sizeof(BAND)); mem_numvar(TYPE_BAND,1); } lb = A->lb = min(n-1,lb); ub = A->ub = min(n-1,ub); A->mat = m_get(lb+ub+1,n); return A; } int bd_free(A) BAND *A; { if ( A == (BAND *)NULL || A->lb < 0 || A->ub < 0 ) /* don't trust it */ return (-1); if (A->mat) m_free(A->mat); if (mem_info_is_on()) { mem_bytes(TYPE_BAND,sizeof(BAND),0); mem_numvar(TYPE_BAND,-1); } free((char *)A); return 0; } /* resize band matrix */ BAND *bd_resize(A,new_lb,new_ub,new_n) BAND *A; int new_lb,new_ub,new_n; { int lb,ub,i,j,l,shift,umin; Real **Av; if (new_lb < 0 || new_ub < 0 || new_n <= 0) error(E_NEG,"bd_resize"); if ( ! A ) return bd_get(new_lb,new_ub,new_n); if ( A->lb+A->ub+1 > A->mat->m ) error(E_INTERN,"bd_resize"); if ( A->lb == new_lb && A->ub == new_ub && A->mat->n == new_n ) return A; lb = A->lb; ub = A->ub; Av = A->mat->me; umin = min(ub,new_ub); /* ensure that unused triangles at edges are zero'd */ for ( i = 0; i < lb; i++ ) for ( j = A->mat->n - lb + i; j < A->mat->n; j++ ) Av[i][j] = 0.0; for ( i = lb+1,l=1; l <= umin; i++,l++ ) for ( j = 0; j < l; j++ ) Av[i][j] = 0.0; new_lb = A->lb = min(new_lb,new_n-1); new_ub = A->ub = min(new_ub,new_n-1); A->mat = m_resize(A->mat,new_lb+new_ub+1,new_n); Av = A->mat->me; /* if new_lb != lb then move the rows to get the main diag in the new_lb row */ if (new_lb > lb) { shift = new_lb-lb; for (i=lb+umin, l=i+shift; i >= 0; i--,l--) MEM_COPY(Av[i],Av[l],new_n*sizeof(Real)); for (l=shift-1; l >= 0; l--) __zero__(Av[l],new_n); } else if (new_lb < lb) { shift = lb - new_lb; for (i=shift, l=0; i <= lb+umin; i++,l++) MEM_COPY(Av[i],Av[l],new_n*sizeof(Real)); for (i=lb+umin+1; i <= new_lb+new_ub; i++) __zero__(Av[i],new_n); } return A; } BAND *bd_copy(A,B) BAND *A,*B; { int lb,ub,i,j,n; if ( !A ) error(E_NULL,"bd_copy"); if (A == B) return B; n = A->mat->n; if ( !B ) B = bd_get(A->lb,A->ub,n); else if (B->lb != A->lb || B->ub != A->ub || B->mat->n != n ) B = bd_resize(B,A->lb,A->ub,n); if (A->mat == B->mat) return B; ub = B->ub = A->ub; lb = B->lb = A->lb; for ( i=0, j=n-lb; i <= lb; i++, j++ ) MEM_COPY(A->mat->me[i],B->mat->me[i],j*sizeof(Real)); for ( i=lb+1, j=1; i <= lb+ub; i++, j++ ) MEM_COPY(A->mat->me[i]+j,B->mat->me[i]+j,(n - j)*sizeof(Real)); return B; } /* copy band matrix to a square matrix */ MAT *band2mat(bA,A) BAND *bA; MAT *A; { int i,j,l,n,n1; int lb, ub; Real **bmat; if ( !bA || !A) error(E_NULL,"band2mat"); if ( bA->mat == A ) error(E_INSITU,"band2mat"); ub = bA->ub; lb = bA->lb; n = bA->mat->n; n1 = n-1; bmat = bA->mat->me; A = m_resize(A,n,n); m_zero(A); for (j=0; j < n; j++) for (i=min(n1,j+lb),l=lb+j-i; i >= max(0,j-ub); i--,l++) A->me[i][j] = bmat[l][j]; return A; } /* copy a square matrix to a band matrix with lb subdiagonals and ub superdiagonals */ BAND *mat2band(A,lb,ub,bA) BAND *bA; MAT *A; int lb, ub; { int i, j, l, n1; Real **bmat; if (! A || ! bA) error(E_NULL,"mat2band"); if (ub < 0 || lb < 0) error(E_SIZES,"mat2band"); if (bA->mat == A) error(E_INSITU,"mat2band"); n1 = A->n-1; lb = min(n1,lb); ub = min(n1,ub); bA = bd_resize(bA,lb,ub,n1+1); bmat = bA->mat->me; for (j=0; j <= n1; j++) for (i=min(n1,j+lb),l=lb+j-i; i >= max(0,j-ub); i--,l++) bmat[l][j] = A->me[i][j]; return bA; } /* transposition of matrix in; out - matrix after transposition; can be done in situ */ BAND *bd_transp(in,out) BAND *in, *out; { int i, j, jj, l, k, lb, ub, lub, n, n1; int in_situ; Real **in_v, **out_v; if ( in == (BAND *)NULL || in->mat == (MAT *)NULL ) error(E_NULL,"bd_transp"); lb = in->lb; ub = in->ub; lub = lb+ub; n = in->mat->n; n1 = n-1; in_situ = ( in == out ); if ( ! in_situ ) out = bd_resize(out,ub,lb,n); else { /* only need to swap lb and ub fields */ out->lb = ub; out->ub = lb; } in_v = in->mat->me; if (! in_situ) { int sh_in,sh_out; out_v = out->mat->me; for (i=0, l=lub, k=lb-i; i <= lub; i++,l--,k--) { sh_in = max(-k,0); sh_out = max(k,0); MEM_COPY(&(in_v[i][sh_in]),&(out_v[l][sh_out]), (n-sh_in-sh_out)*sizeof(Real)); /********************************** for (j=n1-sh_out, jj=n1-sh_in; j >= sh_in; j--,jj--) { out_v[l][jj] = in_v[i][j]; } **********************************/ } } else if (ub == lb) { Real tmp; for (i=0, l=lub, k=lb-i; i < lb; i++,l--,k--) { for (j=n1-k, jj=n1; j >= 0; j--,jj--) { tmp = in_v[l][jj]; in_v[l][jj] = in_v[i][j]; in_v[i][j] = tmp; } } } else if (ub > lb) { /* hence i-ub <= 0 & l-lb >= 0 */ int p,pp,lbi; for (i=0, l=lub; i < (lub+1)/2; i++,l--) { lbi = lb-i; for (j=l-lb, jj=0, p=max(-lbi,0), pp = max(l-ub,0); j <= n1; j++,jj++,p++,pp++) { in_v[l][pp] = in_v[i][p]; in_v[i][jj] = in_v[l][j]; } for ( ; p <= n1-max(lbi,0); p++,pp++) in_v[l][pp] = in_v[i][p]; } if (lub%2 == 0) { /* shift only */ i = lub/2; for (j=max(i-lb,0), jj=0; jj <= n1-ub+i; j++,jj++) in_v[i][jj] = in_v[i][j]; } } else { /* ub < lb, hence ub-l <= 0 & lb-i >= 0 */ int p,pp,ubi; for (i=0, l=lub; i < (lub+1)/2; i++,l--) { ubi = i-ub; for (j=n1-max(lb-l,0), jj=n1-max(-ubi,0), p=n1-lb+i, pp=n1; p >= 0; j--, jj--, pp--, p--) { in_v[i][jj] = in_v[l][j]; in_v[l][pp] = in_v[i][p]; } for ( ; jj >= max(ubi,0); j--, jj--) in_v[i][jj] = in_v[l][j]; } if (lub%2 == 0) { /* shift only */ i = lub/2; for (j=n1-lb+i, jj=n1-max(ub-i,0); j >= 0; j--, jj--) in_v[i][jj] = in_v[i][j]; } } return out; } /* bdLUfactor -- gaussian elimination with partial pivoting -- on entry, the matrix A in band storage with elements in rows 0 to lb+ub; The jth column of A is stored in the jth column of band A (bA) as follows: bA->mat->me[lb+j-i][j] = A->me[i][j] for max(0,j-lb) <= i <= min(A->n-1,j+ub); -- on exit: U is stored as an upper triangular matrix with lb+ub superdiagonals in rows lb to 2*lb+ub, and the matrix L is stored in rows 0 to lb-1. Matrix U is permuted, whereas L is not permuted !!! Therefore we save some memory. */ BAND *bdLUfactor(bA,pivot) BAND *bA; PERM *pivot; { int i, i_max; int j, k, l, n, n1, lb, ub, lub, k_end, k_lub; int shift; Real **bA_v; Real max1, temp; if ( bA==(BAND *)NULL || pivot==(PERM *)NULL ) error(E_NULL,"bdLUfactor"); lb = bA->lb; ub = bA->ub; lub = lb+ub; n = bA->mat->n; n1 = n-1; lub = lb+ub; if ( pivot->size != n ) error(E_SIZES,"bdLUfactor"); /* initialise pivot with identity permutation */ for ( i=0; i < n; i++ ) pivot->pe[i] = i; /* extend band matrix */ /* extended part is filled with zeros */ bA = bd_resize(bA,lb,min(n1,lub),n); bA_v = bA->mat->me; /* main loop */ for ( k=0; k < n1; k++ ) { k_end = max(0,lb+k-n1); k_lub = min(k+lub,n1); /* find the best pivot row */ max1 = 0.0; i_max = -1; for ( i=lb; i >= k_end; i-- ) { temp = fabs(bA_v[i][k]); if ( temp > max1 ) { max1 = temp; i_max = i; } } /* if no pivot then ignore column k... */ if ( i_max == -1 ) continue; /* do we pivot ? */ if ( i_max != lb ) /* yes we do... */ { /* save transposition using non-shifted indices */ shift = lb-i_max; px_transp(pivot,k+shift,k); for ( i=lb, j=k; j <= k_lub; i++,j++ ) { temp = bA_v[i][j]; bA_v[i][j] = bA_v[i-shift][j]; bA_v[i-shift][j] = temp; } } /* row operations */ for ( i=lb-1; i >= k_end; i-- ) { temp = bA_v[i][k] /= bA_v[lb][k]; shift = lb-i; for ( j=k+1,l=i+1; j <= k_lub; l++,j++ ) bA_v[l][j] -= temp*bA_v[l+shift][j]; } } return bA; } /* bdLUsolve -- given an LU factorisation in bA, solve bA*x=b */ /* pivot is changed upon return */ VEC *bdLUsolve(bA,pivot,b,x) BAND *bA; PERM *pivot; VEC *b,*x; { int i,j,l,n,n1,pi,lb,ub,jmin, maxj; Real c; Real **bA_v; if ( bA==(BAND *)NULL || b==(VEC *)NULL || pivot==(PERM *)NULL ) error(E_NULL,"bdLUsolve"); if ( bA->mat->n != b->dim || bA->mat->n != pivot->size) error(E_SIZES,"bdLUsolve"); lb = bA->lb; ub = bA->ub; n = b->dim; n1 = n-1; bA_v = bA->mat->me; x = v_resize(x,b->dim); px_vec(pivot,b,x); /* solve Lx = b; implicit diagonal = 1 L is not permuted, therefore it must be permuted now */ px_inv(pivot,pivot); for (j=0; j < n; j++) { jmin = j+1; c = x->ve[j]; maxj = max(0,j+lb-n1); for (i=jmin,l=lb-1; l >= maxj; i++,l--) { if ( (pi = pivot->pe[i]) < jmin) pi = pivot->pe[i] = pivot->pe[pi]; x->ve[pi] -= bA_v[l][j]*c; } } /* solve Ux = b; explicit diagonal */ x->ve[n1] /= bA_v[lb][n1]; for (i=n-2; i >= 0; i--) { c = x->ve[i]; for (j=min(n1,i+ub), l=lb+j-i; j > i; j--,l--) c -= bA_v[l][j]*x->ve[j]; x->ve[i] = c/bA_v[lb][i]; } return (x); } /* LDLfactor -- L.D.L' factorisation of A in-situ; A is a band matrix it works using only lower bandwidth & main diagonal so it is possible to set A->ub = 0 */ BAND *bdLDLfactor(A) BAND *A; { int i,j,k,n,n1,lb,ki,jk,ji,lbkm,lbkp; Real **Av; Real c, cc; if ( ! A ) error(E_NULL,"bdLDLfactor"); if (A->lb == 0) return A; lb = A->lb; n = A->mat->n; n1 = n-1; Av = A->mat->me; for (k=0; k < n; k++) { lbkm = lb-k; lbkp = lb+k; /* matrix D */ c = Av[lb][k]; for (j=max(0,-lbkm), jk=lbkm+j; j < k; j++, jk++) { cc = Av[jk][j]; c -= Av[lb][j]*cc*cc; } if (c == 0.0) error(E_SING,"bdLDLfactor"); Av[lb][k] = c; /* matrix L */ for (i=min(n1,lbkp), ki=lbkp-i; i > k; i--,ki++) { c = Av[ki][k]; for (j=max(0,i-lb), ji=lb+j-i, jk=lbkm+j; j < k; j++, ji++, jk++) c -= Av[lb][j]*Av[ji][j]*Av[jk][j]; Av[ki][k] = c/Av[lb][k]; } } return A; } /* solve A*x = b, where A is factorized by Choleski LDL^T factorization */ VEC *bdLDLsolve(A,b,x) BAND *A; VEC *b, *x; { int i,j,l,n,n1,lb,ilb; Real **Av, *Avlb; Real c; if ( ! A || ! b ) error(E_NULL,"bdLDLsolve"); if ( A->mat->n != b->dim ) error(E_SIZES,"bdLDLsolve"); n = A->mat->n; n1 = n-1; x = v_resize(x,n); lb = A->lb; Av = A->mat->me; Avlb = Av[lb]; /* solve L*y = b */ x->ve[0] = b->ve[0]; for (i=1; i < n; i++) { ilb = i-lb; c = b->ve[i]; for (j=max(0,ilb), l=j-ilb; j < i; j++,l++) c -= Av[l][j]*x->ve[j]; x->ve[i] = c; } /* solve D*z = y */ for (i=0; i < n; i++) x->ve[i] /= Avlb[i]; /* solve L^T*x = z */ for (i=n-2; i >= 0; i--) { ilb = i+lb; c = x->ve[i]; for (j=min(n1,ilb), l=ilb-j; j > i; j--,l++) c -= Av[l][i]*x->ve[j]; x->ve[i] = c; } return x; } /* ****************************************************** This function is a contribution from Ruediger Franke. His e-mail addres is: Ruediger.Franke@rz.tu-ilmenau.de ****************************************************** */ /* bd_mv_mlt -- * computes out = A * x * may not work in situ (x != out) */ VEC *bd_mv_mlt(BAND *A, VEC *x, VEC *out) { int i, j, j_end, k; int start_idx, end_idx; int n, m, lb, ub; Real **A_me; Real *x_ve; Real sum; if (!A || !x) error(E_NULL,"bd_mv_mlt"); if (x->dim != A->mat->n) error(E_SIZES,"bd_mv_mlt"); if (!out || out->dim != A->mat->n) out = v_resize(out, A->mat->n); if (out == x) error(E_INSITU,"bd_mv_mlt"); n = A->mat->n; m = A->mat->m; lb = A->lb; ub = A->ub; A_me = A->mat->me; start_idx = lb; end_idx = m + n-1 - ub; for (i=0; ive + k; sum = 0.0; for (; j < j_end; j++, k++) sum += A_me[j][k] * *x_ve++; out->ve[i] = sum; } return out; } /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* pxop.c 1.5 12/03/87 */ /********************************************************************** Note: A permutation is often interpreted as a matrix (i.e. a permutation matrix). A permutation px represents a permutation matrix P where P[i][j] == 1 if and only if px->pe[i] == j **********************************************************************/ /* px_inv -- invert permutation -- in situ -- taken from ACM Collected Algorithms #250 */ PERM *px_inv(px,out) PERM *px, *out; { int i, j, k, n, *p; out = px_copy(px, out); n = out->size; p = (int *)(out->pe); for ( n--; n>=0; n-- ) { i = p[n]; if ( i < 0 ) p[n] = -1 - i; else if ( i != n ) { k = n; while (TRUE) { if ( i < 0 || i >= out->size ) error(E_BOUNDS,"px_inv"); j = p[i]; p[i] = -1 - k; if ( j == n ) { p[n] = i; break; } k = i; i = j; } } } return out; } /* px_mlt -- permutation multiplication (composition) */ PERM *px_mlt(px1,px2,out) PERM *px1,*px2,*out; { int i,size; if ( px1==(PERM *)NULL || px2==(PERM *)NULL ) error(E_NULL,"px_mlt"); if ( px1->size != px2->size ) error(E_SIZES,"px_mlt"); if ( px1 == out || px2 == out ) error(E_INSITU,"px_mlt"); if ( out==(PERM *)NULL || out->size < px1->size ) out = px_resize(out,px1->size); size = px1->size; for ( i=0; ipe[i] >= size ) error(E_BOUNDS,"px_mlt"); else out->pe[i] = px1->pe[px2->pe[i]]; return out; } /* px_vec -- permute vector */ VEC *px_vec(px,vector,out) PERM *px; VEC *vector,*out; { int old_i, i, size, start; Real tmp; if ( px==(PERM *)NULL || vector==(VEC *)NULL ) error(E_NULL,"px_vec"); if ( px->size > vector->dim ) error(E_SIZES,"px_vec"); if ( out==(VEC *)NULL || out->dim < vector->dim ) out = v_resize(out,vector->dim); size = px->size; if ( size == 0 ) return v_copy(vector,out); if ( out != vector ) { for ( i=0; ipe[i] >= size ) error(E_BOUNDS,"px_vec"); else out->ve[i] = vector->ve[px->pe[i]]; } else { /* in situ algorithm */ start = 0; while ( start < size ) { old_i = start; i = px->pe[old_i]; if ( i >= size ) { start++; continue; } tmp = vector->ve[start]; while ( TRUE ) { vector->ve[old_i] = vector->ve[i]; px->pe[old_i] = i+size; old_i = i; i = px->pe[old_i]; if ( i >= size ) break; if ( i == start ) { vector->ve[old_i] = tmp; px->pe[old_i] = i+size; break; } } start++; } for ( i = 0; i < size; i++ ) if ( px->pe[i] < size ) error(E_BOUNDS,"px_vec"); else px->pe[i] = px->pe[i]-size; } return out; } /* pxinv_vec -- apply the inverse of px to x, returning the result in out */ VEC *pxinv_vec(px,x,out) PERM *px; VEC *x, *out; { int i, size; if ( ! px || ! x ) error(E_NULL,"pxinv_vec"); if ( px->size > x->dim ) error(E_SIZES,"pxinv_vec"); /* if ( x == out ) error(E_INSITU,"pxinv_vec"); */ if ( ! out || out->dim < x->dim ) out = v_resize(out,x->dim); size = px->size; if ( size == 0 ) return v_copy(x,out); if ( out != x ) { for ( i=0; ipe[i] >= size ) error(E_BOUNDS,"pxinv_vec"); else out->ve[px->pe[i]] = x->ve[i]; } else { /* in situ algorithm --- cheat's way out */ px_inv(px,px); px_vec(px,x,out); px_inv(px,px); } return out; } /* px_transp -- transpose elements of permutation -- Really multiplying a permutation by a transposition */ PERM *px_transp(px,i1,i2) PERM *px; /* permutation to transpose */ int i1,i2; /* elements to transpose */ { int temp; if ( px==(PERM *)NULL ) error(E_NULL,"px_transp"); if ( i1 < px->size && i2 < px->size ) { temp = px->pe[i1]; px->pe[i1] = px->pe[i2]; px->pe[i2] = temp; } return px; } /* myqsort -- a cheap implementation of Quicksort on integers -- returns number of swaps */ static int myqsort(int * a, int num) { int i, j, tmp, v; int numswaps; numswaps = 0; if ( num <= 1 ) return 0; i = 0; j = num; v = a[0]; for ( ; ; ) { while ( a[++i] < v ) ; while ( a[--j] > v ) ; if ( i >= j ) break; tmp = a[i]; a[i] = a[j]; a[j] = tmp; numswaps++; } tmp = a[0]; a[0] = a[j]; a[j] = tmp; if ( j != 0 ) numswaps++; numswaps += myqsort(&a[0],j); numswaps += myqsort(&a[j+1],num-(j+1)); return numswaps; } /* px_sign -- compute the ``sign'' of a permutation = +/-1 where px is the product of an even/odd # transpositions */ int px_sign(px) PERM *px; { int numtransp; PERM *px2; if ( px==(PERM *)NULL ) error(E_NULL,"px_sign"); px2 = px_copy(px,PNULL); numtransp = myqsort(px2->pe,(int)px2->size); px_free(px2); return ( numtransp % 2 ) ? -1 : 1; } /* px_cols -- permute columns of matrix A; out = A.px' -- May NOT be in situ */ MAT *px_cols(px,A,out) PERM *px; MAT *A, *out; { int i, j, m, n, px_j; Real **A_me, **out_me; #ifdef ANSI_C MAT *m_get(int, int); #else extern MAT *m_get(); #endif if ( ! A || ! px ) error(E_NULL,"px_cols"); if ( px->size != A->n ) error(E_SIZES,"px_cols"); if ( A == out ) error(E_INSITU,"px_cols"); m = A->m; n = A->n; if ( ! out || out->m != m || out->n != n ) out = m_get(m,n); A_me = A->me; out_me = out->me; for ( j = 0; j < n; j++ ) { px_j = px->pe[j]; if ( px_j >= n ) error(E_BOUNDS,"px_cols"); for ( i = 0; i < m; i++ ) out_me[i][px_j] = A_me[i][j]; } return out; } /* px_rows -- permute columns of matrix A; out = px.A -- May NOT be in situ */ MAT *px_rows(px,A,out) PERM *px; MAT *A, *out; { int i, j, m, n, px_i; Real **A_me, **out_me; #ifdef ANSI_C MAT *m_get(int, int); #else extern MAT *m_get(); #endif if ( ! A || ! px ) error(E_NULL,"px_rows"); if ( px->size != A->m ) error(E_SIZES,"px_rows"); if ( A == out ) error(E_INSITU,"px_rows"); m = A->m; n = A->n; if ( ! out || out->m != m || out->n != n ) out = m_get(m,n); A_me = A->me; out_me = out->me; for ( i = 0; i < m; i++ ) { px_i = px->pe[i]; if ( px_i >= m ) error(E_BOUNDS,"px_rows"); for ( j = 0; j < n; j++ ) out_me[i][j] = A_me[px_i][j]; } return out; } /************************************************************************** ** ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved. ** ** Meschach Library ** ** This Meschach Library is provided "as is" without any express ** or implied warranty of any kind with respect to this software. ** In particular the authors shall not be liable for any direct, ** indirect, special, incidental or consequential damages arising ** in any way from use of the software. ** ** Everyone is granted permission to copy, modify and redistribute this ** Meschach Library, provided: ** 1. All copies contain this copyright notice. ** 2. All modified copies shall carry a notice stating who ** made the last modification and the date of such modification. ** 3. No charge is made for this software or works derived from it. ** This clause shall not be construed as constraining other software ** distributed on the same medium as this software, nor is a ** distribution fee considered a charge. ** ***************************************************************************/ /* _m_copy -- copies matrix into new area */ MAT *_m_copy(in,out,i0,j0) MAT *in,*out; int i0,j0; { int i /* ,j */; if ( in==MNULL ) error(E_NULL,"_m_copy"); if ( in==out ) return (out); if ( out==MNULL || out->m < in->m || out->n < in->n ) out = m_resize(out,in->m,in->n); for ( i=i0; i < in->m; i++ ) MEM_COPY(&(in->me[i][j0]),&(out->me[i][j0]), (in->n - j0)*sizeof(Real)); /* for ( j=j0; j < in->n; j++ ) out->me[i][j] = in->me[i][j]; */ return (out); } /* _v_copy -- copies vector into new area */ VEC *_v_copy(in,out,i0) VEC *in,*out; int i0; { /* int i,j; */ if ( in==VNULL ) error(E_NULL,"_v_copy"); if ( i0 < 0 || in->dim < i0 ) error(E_RANGE,"_v_copy"); if ( in==out ) return (out); if ( out==VNULL || out->dim < in->dim ) out = v_resize(out,in->dim); MEM_COPY(&(in->ve[i0]),&(out->ve[i0]),(in->dim - i0)*sizeof(Real)); /* for ( i=i0; i < in->dim; i++ ) out->ve[i] = in->ve[i]; */ return (out); } /* px_copy -- copies permutation 'in' to 'out' */ PERM *px_copy(in,out) PERM *in,*out; { /* int i; */ if ( in == PNULL ) error(E_NULL,"px_copy"); if ( in == out ) return out; if ( out == PNULL || out->size != in->size ) out = px_resize(out,in->size); MEM_COPY(in->pe,out->pe,in->size*sizeof(int)); /* for ( i = 0; i < in->size; i++ ) out->pe[i] = in->pe[i]; */ return out; } /* The .._move() routines are for moving blocks of memory around within Meschach data structures and for re-arranging matrices, vectors etc. */ /* m_move -- copies selected pieces of a matrix -- moves the m0 x n0 submatrix with top-left cor-ordinates (i0,j0) to the corresponding submatrix of out with top-left co-ordinates (i1,j1) -- out is resized (& created) if necessary */ MAT *m_move(in,i0,j0,m0,n0,out,i1,j1) MAT *in, *out; int i0, j0, m0, n0, i1, j1; { int i; if ( ! in ) error(E_NULL,"m_move"); if ( i0 < 0 || j0 < 0 || i1 < 0 || j1 < 0 || m0 < 0 || n0 < 0 || i0+m0 > in->m || j0+n0 > in->n ) error(E_BOUNDS,"m_move"); if ( ! out ) out = m_resize(out,i1+m0,j1+n0); else if ( i1+m0 > out->m || j1+n0 > out->n ) out = m_resize(out,max(out->m,i1+m0),max(out->n,j1+n0)); for ( i = 0; i < m0; i++ ) MEM_COPY(&(in->me[i0+i][j0]),&(out->me[i1+i][j1]), n0*sizeof(Real)); return out; } /* v_move -- copies selected pieces of a vector -- moves the length dim0 subvector with initial index i0 to the corresponding subvector of out with initial index i1 -- out is resized if necessary */ VEC *v_move(in,i0,dim0,out,i1) VEC *in, *out; int i0, dim0, i1; { if ( ! in ) error(E_NULL,"v_move"); if ( i0 < 0 || dim0 < 0 || i1 < 0 || i0+dim0 > in->dim ) error(E_BOUNDS,"v_move"); if ( (! out) || i1+dim0 > out->dim ) out = v_resize(out,i1+dim0); MEM_COPY(&(in->ve[i0]),&(out->ve[i1]),dim0*sizeof(Real)); return out; } /* mv_move -- copies selected piece of matrix to a vector -- moves the m0 x n0 submatrix with top-left co-ordinate (i0,j0) to the subvector with initial index i1 (and length m0*n0) -- rows are copied contiguously -- out is resized if necessary */ VEC *mv_move(in,i0,j0,m0,n0,out,i1) MAT *in; VEC *out; int i0, j0, m0, n0, i1; { int dim1, i; if ( ! in ) error(E_NULL,"mv_move"); if ( i0 < 0 || j0 < 0 || m0 < 0 || n0 < 0 || i1 < 0 || i0+m0 > in->m || j0+n0 > in->n ) error(E_BOUNDS,"mv_move"); dim1 = m0*n0; if ( (! out) || i1+dim1 > out->dim ) out = v_resize(out,i1+dim1); for ( i = 0; i < m0; i++ ) MEM_COPY(&(in->me[i0+i][j0]),&(out->ve[i1+i*n0]),n0*sizeof(Real)); return out; } /* vm_move -- copies selected piece of vector to a matrix -- moves the subvector with initial index i0 and length m1*n1 to the m1 x n1 submatrix with top-left co-ordinate (i1,j1) -- copying is done by rows -- out is resized if necessary */ MAT *vm_move(in,i0,out,i1,j1,m1,n1) VEC *in; MAT *out; int i0, i1, j1, m1, n1; { int dim0, i; if ( ! in ) error(E_NULL,"vm_move"); if ( i0 < 0 || i1 < 0 || j1 < 0 || m1 < 0 || n1 < 0 || i0+m1*n1 > in->dim ) error(E_BOUNDS,"vm_move"); if ( ! out ) out = m_resize(out,i1+m1,j1+n1); else out = m_resize(out,max(i1+m1,out->m),max(j1+n1,out->n)); dim0 = m1*n1; for ( i = 0; i < m1; i++ ) MEM_COPY(&(in->ve[i0+i*n1]),&(out->me[i1+i][j1]),n1*sizeof(Real)); return out; } /* * getdelim.c -- this function copies the GNU glibc function * int getdelim( char **ppline, size_t* pscLen, int sc, FILE* pF); * which reads a malloc limited string into the buffer * * Pre condition: * *ppline == NULL && *psclen == 0 : call malloc for new memory * *ppline != NULL && *psclen >= 0 : realloc memory * sc = '/n' or some other delimiter character * pf = stdin or some other file pointer * * Post condition: * Returns number of characters read or EOF * if it cannot find a delimiter, it reads until EOF, and returns the * number of characters read, and an EOF on the next read. * *pline = pointer to null delimited string / read buffer * *psclen = new length of read buffer memory * * dave forrest 1998 July 13 */ #define GETDELIM_EXT 128 /* use -DTESTGETDELIM for a test program */ int getdelim( char **ppline, size_t* pscLen, int sc, FILE* pF); int getdelim( char **ppLine, size_t* pscLen, int sc, FILE* pF){ unsigned int scRead; int scX; char* pBuff; scRead = 0; if(ppLine ==NULL || pscLen == NULL) return -1; /* bad setup */ if(*ppLine ==NULL && *pscLen != 0 ) return -1; /* bad call */ if(pF == NULL) return -1; /* bad file pointer */ while ( (scX = fgetc(pF)) != EOF) { if (scRead+2 > *pscLen) { /* Extention required */ if ( (pBuff = (char *)realloc(*ppLine,*pscLen+GETDELIM_EXT)) == NULL ) return -1; else { *ppLine = pBuff ; *pscLen+=GETDELIM_EXT; } } (*ppLine)[scRead] = scX; scRead++; if(scX == sc) break; /* found the delimiter */ } (*ppLine)[scRead] = '\0'; if(scRead == 0) return -1; /* End of File */ return scRead ; } #ifdef TESTGETDELIM int main(){ char* plinebuffer ; size_t linebuffersize; int length; printf("\n\ getdelim.c -- reads a file for a delimited string of arbitrary length \n\ and returns the number of characters read, a pointer to the \n\ malloc'd buffer, and a length of the buffer -- EOF on error \n\ dforrest@virginia.edu 1998 July 13 \n\ \n\ "); printf("test--&pbuff &bufsiz 3 tests for null calls\n"); if( getdelim( &plinebuffer , NULL ,'\n', stdin) == -1) printf("OK -- xxxx NULL\n"); if( getdelim( NULL , &linebuffersize ,'\n', stdin) == -1) printf("OK -- NULL xxxx\n"); linebuffersize = 5; if( getdelim( &plinebuffer, &linebuffersize ,'\n', stdin) == -1) printf("OK -- xxxxx xxxx>5\n"); else printf("bad trapping of null pointer and non-zero buffer length\n"); linebuffersize =0; while ((length = getdelim(&plinebuffer,&linebuffersize,'\n', stdin)) > 0) { printf ("%d :'%s'\n",length, plinebuffer); } printf ("%d : *** returned from getdelim (0 = no more, -1 = EOF)\n",length); return 0; } #endif /* * readfloat.c -- read a string full of floats into an array * */ /* compile with -DTEST to get a test program */ int readfloats(char *pst, float ** ppvcX, int * pscDimMax) { /* pst is the pointer to the string pvcX is a pointer to an array of floats scDimMax is the maximum size of the array */ int sc; /* counter */ int scAte; /* characters consumed by read */ double scX; /* current value */ char * pstRest; /* remainder of string */ void * p; /* pointer */ pstRest = pst; sc = 0; /* count */ if (pst ==NULL || pst[0] == 0) return 0; /* handle null input */ if (*ppvcX == NULL && *pscDimMax > 0 ) return -1 ; /* bad setup */ while ( sscanf(pstRest,"%lf%n",&scX,&scAte) >=1 ) { sc++; if (sc >*pscDimMax) { /* array not big enough */ if( (p=realloc( (*ppvcX), (*pscDimMax+20)*sizeof(float)))==NULL){ perror("realloc failed"); exit(1); }else{ *ppvcX=p; *pscDimMax += 20; #if DEBUG > 5 printf("Extended array readvec read %d values\nfrom %s\n leaving %s",sc,pst,pstRest); #endif } } (*ppvcX)[sc-1] = (float)scX; pstRest+=scAte; } #if DEBUG > 5 printf("readvec read %d values\nfrom:%s\nleaving:%s",sc,pst,pstRest); #endif return (sc) ; } /* +++Date last modified: 05-Jul-1997 */ /* ======================================================================= CFG.c Generic configuration file handler. A. Reitsma, Delft, The Netherlands. v1.00 94-07-09 Public Domain. David Forrest, Charlottesville, USA drf5n@virginia.edu v1.1 98-08-01 + Accepts multi-line quotes using "'` + Ignores # and ; as comments + malloc-(un)limited input using GNU's getdelim + changed the return values - = bad + = # of assignments + added case-sensitivity Bugs - Multiple assignments will free() old ones. -variables may be free()'d if clobbered, so don't set ---------------------------------------------------------------------- */ /* #define TEST */ /* #define DEBUG */ enum RetVal { ERR_FOPEN=-1, ERR_MEM=-2, ERR_QUOTE=-3, ERR_FCLOSE=-4 }; /* end of enum */ int CfgRead( char * Filename, struct CfgStrings * CfgInfo ){ char *psBuffer = NULL; size_t cBufferMax = 0; int chQuote; int cValues = 0; char *psQuoteBuffer = NULL; size_t cQuoteMax = 0; char * pstWork ; char * psStart ; char * CfgName ; char * CfgData ; char * psQuote ; int cQuote, cCfgData; short fInQuote = 0; struct CfgStrings * Cfg ; FILE * CfgFile ; CfgFile = fopen( Filename, "r" ); if( NULL == CfgFile ) { return ERR_FOPEN ; } while( getdelim( &psBuffer, &cBufferMax, '\n', CfgFile ) > 0) { /* While reading a \n terminated line... */ DebugCode( printf (" getdelim:%s",psBuffer); ) /* clip off Leading Whitespace & dump the Comments/Blanks */ psStart = psBuffer; while( isspace( *psStart )) psStart++; /* if( (index("#;\0\r\n" ,*psStart )!=0) || 0 == strlen( psStart ))*/ if( (strchr("#;\0\r\n" ,*psStart )!=0) || 0 == strlen( psStart )) continue; /* Eat the empty and comment lines */ DebugCode( printf ("Using:%s",psStart); ) CfgName = strtok( psStart, " =" ); if( NULL != CfgName ) {/* We are Reading a line with a name */ fInQuote = 0; CfgData = strtok( NULL, "" ); DebugCode(printf ("CfgName=:%s:",CfgName); ) while( isspace( *CfgData )) CfgData++; if( '=' == *CfgData ) CfgData++; while( isspace( *CfgData )) CfgData++; /* if ( index("`'\"",CfgData[0] ) != NULL) */ if ( strchr("`'\"",CfgData[0] ) != NULL) { /* Got a quote character */ chQuote = CfgData[0]; psQuote= ++CfgData; DebugCode(printf ("quote:%c:",chQuote);) if ( (pstWork=strchr(psQuote,chQuote)) != NULL) { *pstWork = '\0'; } else{ /* get up until the next quote */ cQuote = getdelim(&psQuoteBuffer, &cQuoteMax, chQuote, CfgFile); if (cQuote < 1) { if (fclose( CfgFile ) != 0 ){return ERR_FCLOSE;} return ERR_QUOTE; } if(psQuoteBuffer[cQuote-1] !=chQuote) return ERR_QUOTE; psQuoteBuffer[cQuote-1] = '\0'; fInQuote = 1; } } /* Quotes are handled */ else { /* strip off trailing comments */ if ((pstWork = strchr(CfgData,';')) != NULL) { *(pstWork--)= '\0'; while(isspace(*pstWork)) *(pstWork--)='\0'; } if ((pstWork = strchr(CfgData,'#')) != NULL) { *(pstWork--)= '\0'; while(isspace(*pstWork)) *(pstWork--)='\0'; while((*pstWork) == '\r') *(pstWork--)='\0'; } } cCfgData=strlen(CfgData)+1; DebugCode(printf ("Data:%s:\n",CfgData); if(fInQuote) printf("--and:%s:\n",psQuoteBuffer); ) /* look for matching 'name' */ Cfg = CfgInfo ; while( (NULL != Cfg->name) && (0 != strcmp( Cfg->name, CfgName ))) Cfg++; /* duplicate the data if the name is found. */ if( NULL != Cfg->name ) { /* Our name matches a name in the list */ cValues++; free(Cfg->data); /* Cfg->data = strdup( CfgData ); */ Cfg->data = malloc((unsigned)cCfgData); if( NULL == Cfg->data ) { /* strdup (malloc) failed */ fclose( CfgFile ); return ERR_MEM; } else { memcpy(Cfg->data,CfgData,(unsigned)cCfgData); } if (fInQuote) { /* More stuff in the Quoted buffer */ Cfg->data = (char *) realloc(Cfg->data, strlen(Cfg->data)+cQuote); if( NULL == Cfg->data ) { /* realloc failed */ fclose( CfgFile ); return ERR_MEM; } strcat(Cfg->data,psQuoteBuffer); } } } } free(psBuffer); free(psQuoteBuffer); if (fclose( CfgFile ) != 0 ){ return ERR_FCLOSE;}; return cValues ; } /* End Of cfg */ /* ==== CFG.c end ===================================================== */ /* getopt.c from the GNU people */ /* Getopt for GNU. NOTE: getopt is now part of the C library, so if you don't know what "Keep this file name-space clean" means, talk to roland@gnu.ai.mit.edu before changing it! Copyright (C) 1987, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97 Free Software Foundation, Inc. The GNU C Library is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. The GNU C Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more details. You should have received a copy of the GNU Library General Public License along with the GNU C Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ /* This tells Alpha OSF/1 not to define a getopt prototype in . Ditto for AIX 3.2 and . */ #ifndef _NO_PROTO #define _NO_PROTO #endif #ifdef HAVE_CONFIG_H #include #endif #if !defined (__STDC__) || !__STDC__ /* This is a separate conditional since some stdc systems reject `defined (const)'. */ #ifndef const #define const #endif #endif /* Comment out all this code if we are using the GNU C Library, and are not actually compiling the library itself. This code is part of the GNU C Library, but also included in many other GNU distributions. Compiling and linking in this code is a waste when using the GNU C library (especially if it is a shared library). Rather than having every GNU program understand `configure --with-gnu-libc' and omit the object files, it is simpler to just do this in the source for each such file. */ #define GETOPT_INTERFACE_VERSION 2 #if !defined (_LIBC) && defined (__GLIBC__) && __GLIBC__ >= 2 #include #if _GNU_GETOPT_INTERFACE_VERSION == GETOPT_INTERFACE_VERSION #define ELIDE_CODE #endif #endif #ifndef ELIDE_CODE /* This needs to come after some library #include to get __GNU_LIBRARY__ defined. */ #ifdef __GNU_LIBRARY__ /* Don't include stdlib.h for non-GNU C libraries because some of them contain conflicting prototypes for getopt. */ #include /*#include */ #endif /* GNU C library. */ #ifdef VMS #include #if HAVE_STRING_H - 0 #include #endif #endif #if defined (WIN32) && !defined (__CYGWIN32__) /* It's not Unix, really. See? Capital letters. */ /* #include */ #define getpid() GetCurrentProcessId() #endif #ifndef _ /* This is for other GNU distributions with internationalized messages. When compiling libc, the _ macro is predefined. */ #ifdef HAVE_LIBINTL_H # include # define _(msgid) gettext (msgid) #else # define _(msgid) (msgid) #endif #endif /* This version of `getopt' appears to the caller like standard Unix `getopt' but it behaves differently for the user, since it allows the user to intersperse the options with the other arguments. As `getopt' works, it permutes the elements of ARGV so that, when it is done, all the options precede everything else. Thus all application programs are extended to handle flexible argument order. Setting the environment variable POSIXLY_CORRECT disables permutation. Then the behavior is completely standard. GNU application programs can use a third alternative mode in which they can distinguish the relative order of options and other arguments. */ /* #include "getopt.h" */ /* For communication from `getopt' to the caller. When `getopt' finds an option that takes an argument, the argument value is returned here. Also, when `ordering' is RETURN_IN_ORDER, each non-option ARGV-element is returned here. */ char *optarg = NULL; /* Index in ARGV of the next element to be scanned. This is used for communication to and from the caller and for communication between successive calls to `getopt'. On entry to `getopt', zero means this is the first call; initialize. When `getopt' returns -1, this is the index of the first of the non-option elements that the caller should itself scan. Otherwise, `optind' communicates from one call to the next how much of ARGV has been scanned so far. */ /* 1003.2 says this must be 1 before any call. */ int optind = 1; /* Formerly, initialization of getopt depended on optind==0, which causes problems with re-calling getopt as programs generally don't know that. */ int __getopt_initialized = 0; /* The next char to be scanned in the option-element in which the last option character we returned was found. This allows us to pick up the scan where we left off. If this is zero, or a null string, it means resume the scan by advancing to the next ARGV-element. */ static char *nextchar; /* Callers store zero here to inhibit the error message for unrecognized options. */ int opterr = 1; /* Set to an option character which was unrecognized. This must be initialized on some systems to avoid linking in the system's own getopt implementation. */ int optopt = '?'; /* Describe how to deal with options that follow non-option ARGV-elements. If the caller did not specify anything, the default is REQUIRE_ORDER if the environment variable POSIXLY_CORRECT is defined, PERMUTE otherwise. REQUIRE_ORDER means don't recognize them as options; stop option processing when the first non-option is seen. This is what Unix does. This mode of operation is selected by either setting the environment variable POSIXLY_CORRECT, or using `+' as the first character of the list of option characters. PERMUTE is the default. We permute the contents of ARGV as we scan, so that eventually all the non-options are at the end. This allows options to be given in any order, even with programs that were not written to expect this. RETURN_IN_ORDER is an option available to programs that were written to expect options and other ARGV-elements in any order and that care about the ordering of the two. We describe each non-option ARGV-element as if it were the argument of an option with character code 1. Using `-' as the first character of the list of option characters selects this mode of operation. The special argument `--' forces an end of option-scanning regardless of the value of `ordering'. In the case of RETURN_IN_ORDER, only `--' can cause `getopt' to return -1 with `optind' != ARGC. */ static enum { REQUIRE_ORDER, PERMUTE, RETURN_IN_ORDER } ordering; /* Value of POSIXLY_CORRECT environment variable. */ static char *posixly_correct; #ifdef __GNU_LIBRARY__ /* We want to avoid inclusion of string.h with non-GNU libraries because there are many ways it can cause trouble. On some systems, it contains special magic macros that don't work in GCC. */ #include #define my_index strchr #else /* Avoid depending on library functions or files whose names are inconsistent. */ char *getenv (); static char * my_index (str, chr) const char *str; int chr; { while (*str) { if (*str == chr) return (char *) str; str++; } return 0; } /* If using GCC, we can safely declare strlen this way. If not using GCC, it is ok not to declare it. */ #ifdef __GNUC__ /* Note that Motorola Delta 68k R3V7 comes with GCC but not stddef.h. That was relevant to code that was here before. */ #if !defined (__STDC__) || !__STDC__ /* gcc with -traditional declares the built-in strlen to return int, and has done so at least since version 2.4.5. -- rms. */ extern int strlen (const char *); #endif /* not __STDC__ */ #endif /* __GNUC__ */ #endif /* not __GNU_LIBRARY__ */ /* Handle permutation of arguments. */ /* Describe the part of ARGV that contains non-options that have been skipped. `first_nonopt' is the index in ARGV of the first of them; `last_nonopt' is the index after the last of them. */ static int first_nonopt; static int last_nonopt; #ifdef _LIBC /* Bash 2.0 gives us an environment variable containing flags indicating ARGV elements that should not be considered arguments. */ /* Defined in getopt_init.c */ extern char *__getopt_nonoption_flags; static int nonoption_flags_max_len; static int nonoption_flags_len; static int original_argc; static char *const *original_argv; extern pid_t __libc_pid; /* Make sure the environment variable bash 2.0 puts in the environment is valid for the getopt call we must make sure that the ARGV passed to getopt is that one passed to the process. */ static void __attribute__ ((unused)) store_args_and_env (int argc, char *const *argv) { /* XXX This is no good solution. We should rather copy the args so that we can compare them later. But we must not use malloc(3). */ original_argc = argc; original_argv = argv; } text_set_element (__libc_subinit, store_args_and_env); # define SWAP_FLAGS(ch1, ch2) \ if (nonoption_flags_len > 0) \ { \ char __tmp = __getopt_nonoption_flags[ch1]; \ __getopt_nonoption_flags[ch1] = __getopt_nonoption_flags[ch2]; \ __getopt_nonoption_flags[ch2] = __tmp; \ } #else /* !_LIBC */ # define SWAP_FLAGS(ch1, ch2) #endif /* _LIBC */ /* Exchange two adjacent subsequences of ARGV. One subsequence is elements [first_nonopt,last_nonopt) which contains all the non-options that have been skipped so far. The other is elements [last_nonopt,optind), which contains all the options processed since those non-options were skipped. `first_nonopt' and `last_nonopt' are relocated so that they describe the new indices of the non-options in ARGV after they are moved. */ #if defined (__STDC__) && __STDC__ static void exchange (char **); #endif static void exchange (argv) char **argv; { int bottom = first_nonopt; int middle = last_nonopt; int top = optind; char *tem; /* Exchange the shorter segment with the far end of the longer segment. That puts the shorter segment into the right place. It leaves the longer segment in the right place overall, but it consists of two parts that need to be swapped next. */ #ifdef _LIBC /* First make sure the handling of the `__getopt_nonoption_flags' string can work normally. Our top argument must be in the range of the string. */ if (nonoption_flags_len > 0 && top >= nonoption_flags_max_len) { /* We must extend the array. The user plays games with us and presents new arguments. */ char *new_str = malloc (top + 1); if (new_str == NULL) nonoption_flags_len = nonoption_flags_max_len = 0; else { memcpy (new_str, __getopt_nonoption_flags, nonoption_flags_max_len); memset (&new_str[nonoption_flags_max_len], '\0', top + 1 - nonoption_flags_max_len); nonoption_flags_max_len = top + 1; __getopt_nonoption_flags = new_str; } } #endif while (top > middle && middle > bottom) { if (top - middle > middle - bottom) { /* Bottom segment is the short one. */ int len = middle - bottom; int i; /* Swap it with the top part of the top segment. */ for (i = 0; i < len; i++) { tem = argv[bottom + i]; argv[bottom + i] = argv[top - (middle - bottom) + i]; argv[top - (middle - bottom) + i] = tem; SWAP_FLAGS (bottom + i, top - (middle - bottom) + i); } /* Exclude the moved bottom segment from further swapping. */ top -= len; } else { /* Top segment is the short one. */ int len = top - middle; int i; /* Swap it with the bottom part of the bottom segment. */ for (i = 0; i < len; i++) { tem = argv[bottom + i]; argv[bottom + i] = argv[middle + i]; argv[middle + i] = tem; SWAP_FLAGS (bottom + i, middle + i); } /* Exclude the moved top segment from further swapping. */ bottom += len; } } /* Update records for the slots the non-options now occupy. */ first_nonopt += (optind - last_nonopt); last_nonopt = optind; } /* Initialize the internal data when the first call is made. */ #if defined (__STDC__) && __STDC__ static const char *_getopt_initialize (int, char *const *, const char *); #endif static const char * _getopt_initialize (argc, argv, optstring) int argc; char *const *argv; const char *optstring; { /* Start processing options with ARGV-element 1 (since ARGV-element 0 is the program name); the sequence of previously skipped non-option ARGV-elements is empty. */ first_nonopt = last_nonopt = optind; nextchar = NULL; posixly_correct = getenv ("POSIXLY_CORRECT"); /* Determine how to handle the ordering of options and nonoptions. */ if (optstring[0] == '-') { ordering = RETURN_IN_ORDER; ++optstring; } else if (optstring[0] == '+') { ordering = REQUIRE_ORDER; ++optstring; } else if (posixly_correct != NULL) ordering = REQUIRE_ORDER; else ordering = PERMUTE; #ifdef _LIBC if (posixly_correct == NULL && argc == original_argc && argv == original_argv) { if (nonoption_flags_max_len == 0) { if (__getopt_nonoption_flags == NULL || __getopt_nonoption_flags[0] == '\0') nonoption_flags_max_len = -1; else { const char *orig_str = __getopt_nonoption_flags; int len = nonoption_flags_max_len = strlen (orig_str); if (nonoption_flags_max_len < argc) nonoption_flags_max_len = argc; __getopt_nonoption_flags = (char *) malloc (nonoption_flags_max_len); if (__getopt_nonoption_flags == NULL) nonoption_flags_max_len = -1; else { memcpy (__getopt_nonoption_flags, orig_str, len); memset (&__getopt_nonoption_flags[len], '\0', nonoption_flags_max_len - len); } } } nonoption_flags_len = nonoption_flags_max_len; } else nonoption_flags_len = 0; #endif return optstring; } /* Scan elements of ARGV (whose length is ARGC) for option characters given in OPTSTRING. If an element of ARGV starts with '-', and is not exactly "-" or "--", then it is an option element. The characters of this element (aside from the initial '-') are option characters. If `getopt' is called repeatedly, it returns successively each of the option characters from each of the option elements. If `getopt' finds another option character, it returns that character, updating `optind' and `nextchar' so that the next call to `getopt' can resume the scan with the following option character or ARGV-element. If there are no more option characters, `getopt' returns -1. Then `optind' is the index in ARGV of the first ARGV-element that is not an option. (The ARGV-elements have been permuted so that those that are not options now come last.) OPTSTRING is a string containing the legitimate option characters. If an option character is seen that is not listed in OPTSTRING, return '?' after printing an error message. If you set `opterr' to zero, the error message is suppressed but we still return '?'. If a char in OPTSTRING is followed by a colon, that means it wants an arg, so the following text in the same ARGV-element, or the text of the following ARGV-element, is returned in `optarg'. Two colons mean an option that wants an optional arg; if there is text in the current ARGV-element, it is returned in `optarg', otherwise `optarg' is set to zero. If OPTSTRING starts with `-' or `+', it requests different methods of handling the non-option ARGV-elements. See the comments about RETURN_IN_ORDER and REQUIRE_ORDER, above. Long-named options begin with `--' instead of `-'. Their names may be abbreviated as long as the abbreviation is unique or is an exact match for some defined option. If they have an argument, it follows the option name in the same ARGV-element, separated from the option name by a `=', or else the in next ARGV-element. When `getopt' finds a long-named option, it returns 0 if that option's `flag' field is nonzero, the value of the option's `val' field if the `flag' field is zero. The elements of ARGV aren't really const, because we permute them. But we pretend they're const in the prototype to be compatible with other systems. LONGOPTS is a vector of `struct option' terminated by an element containing a name which is zero. LONGIND returns the index in LONGOPT of the long-named option found. It is only valid when a long-named option has been found by the most recent call. If LONG_ONLY is nonzero, '-' as well as '--' can introduce long-named options. */ int _getopt_internal (argc, argv, optstring, longopts, longind, long_only) int argc; char *const *argv; const char *optstring; const struct option *longopts; int *longind; int long_only; { optarg = NULL; if (optind == 0 || !__getopt_initialized) { if (optind == 0) optind = 1; /* Don't scan ARGV[0], the program name. */ optstring = _getopt_initialize (argc, argv, optstring); __getopt_initialized = 1; } /* Test whether ARGV[optind] points to a non-option argument. Either it does not have option syntax, or there is an environment flag from the shell indicating it is not an option. The later information is only used when the used in the GNU libc. */ #ifdef _LIBC #define NONOPTION_P (argv[optind][0] != '-' || argv[optind][1] == '\0' \ || (optind < nonoption_flags_len \ && __getopt_nonoption_flags[optind] == '1')) #else #define NONOPTION_P (argv[optind][0] != '-' || argv[optind][1] == '\0') #endif if (nextchar == NULL || *nextchar == '\0') { /* Advance to the next ARGV-element. */ /* Give FIRST_NONOPT & LAST_NONOPT rational values if OPTIND has been moved back by the user (who may also have changed the arguments). */ if (last_nonopt > optind) last_nonopt = optind; if (first_nonopt > optind) first_nonopt = optind; if (ordering == PERMUTE) { /* If we have just processed some options following some non-options, exchange them so that the options come first. */ if (first_nonopt != last_nonopt && last_nonopt != optind) exchange ((char **) argv); else if (last_nonopt != optind) first_nonopt = optind; /* Skip any additional non-options and extend the range of non-options previously skipped. */ while (optind < argc && NONOPTION_P) optind++; last_nonopt = optind; } /* The special ARGV-element `--' means premature end of options. Skip it like a null option, then exchange with previous non-options as if it were an option, then skip everything else like a non-option. */ if (optind != argc && !strcmp (argv[optind], "--")) { optind++; if (first_nonopt != last_nonopt && last_nonopt != optind) exchange ((char **) argv); else if (first_nonopt == last_nonopt) first_nonopt = optind; last_nonopt = argc; optind = argc; } /* If we have done all the ARGV-elements, stop the scan and back over any non-options that we skipped and permuted. */ if (optind == argc) { /* Set the next-arg-index to point at the non-options that we previously skipped, so the caller will digest them. */ if (first_nonopt != last_nonopt) optind = first_nonopt; return -1; } /* If we have come to a non-option and did not permute it, either stop the scan or describe it to the caller and pass it by. */ if (NONOPTION_P) { if (ordering == REQUIRE_ORDER) return -1; optarg = argv[optind++]; return 1; } /* We have found another option-ARGV-element. Skip the initial punctuation. */ nextchar = (argv[optind] + 1 + (longopts != NULL && argv[optind][1] == '-')); } /* Decode the current option-ARGV-element. */ /* Check whether the ARGV-element is a long option. If long_only and the ARGV-element has the form "-f", where f is a valid short option, don't consider it an abbreviated form of a long option that starts with f. Otherwise there would be no way to give the -f short option. On the other hand, if there's a long option "fubar" and the ARGV-element is "-fu", do consider that an abbreviation of the long option, just like "--fu", and not "-f" with arg "u". This distinction seems to be the most useful approach. */ if (longopts != NULL && (argv[optind][1] == '-' || (long_only && (argv[optind][2] || !my_index (optstring, argv[optind][1]))))) { char *nameend; const struct option *p; const struct option *pfound = NULL; int exact = 0; int ambig = 0; int indfound = -1; int option_index; for (nameend = nextchar; *nameend && *nameend != '='; nameend++) /* Do nothing. */ ; /* Test all long options for either exact match or abbreviated matches. */ for (p = longopts, option_index = 0; p->name; p++, option_index++) if (!strncmp (p->name, nextchar, (unsigned)(nameend - nextchar))) { if ((unsigned int) (nameend - nextchar) == (unsigned int) strlen (p->name)) { /* Exact match found. */ pfound = p; indfound = option_index; exact = 1; break; } else if (pfound == NULL) { /* First nonexact match found. */ pfound = p; indfound = option_index; } else /* Second or later nonexact match found. */ ambig = 1; } if (ambig && !exact) { if (opterr) fprintf (stderr, _("%s: option `%s' is ambiguous\n"), argv[0], argv[optind]); nextchar += strlen (nextchar); optind++; optopt = 0; return '?'; } if (pfound != NULL) { option_index = indfound; optind++; if (*nameend) { /* Don't test has_arg with >, because some C compilers don't allow it to be used on enums. */ if (pfound->has_arg) optarg = nameend + 1; else { if (opterr) if (argv[optind - 1][1] == '-') /* --option */ fprintf (stderr, _("%s: option `--%s' doesn't allow an argument\n"), argv[0], pfound->name); else /* +option or -option */ fprintf (stderr, _("%s: option `%c%s' doesn't allow an argument\n"), argv[0], argv[optind - 1][0], pfound->name); nextchar += strlen (nextchar); optopt = pfound->val; return '?'; } } else if (pfound->has_arg == 1) { if (optind < argc) optarg = argv[optind++]; else { if (opterr) fprintf (stderr, _("%s: option `%s' requires an argument\n"), argv[0], argv[optind - 1]); nextchar += strlen (nextchar); optopt = pfound->val; return optstring[0] == ':' ? ':' : '?'; } } nextchar += strlen (nextchar); if (longind != NULL) *longind = option_index; if (pfound->flag) { *(pfound->flag) = pfound->val; return 0; } return pfound->val; } /* Can't find it as a long option. If this is not getopt_long_only, or the option starts with '--' or is not a valid short option, then it's an error. Otherwise interpret it as a short option. */ if (!long_only || argv[optind][1] == '-' || my_index (optstring, *nextchar) == NULL) { if (opterr) { if (argv[optind][1] == '-') /* --option */ fprintf (stderr, _("%s: unrecognized option `--%s'\n"), argv[0], nextchar); else /* +option or -option */ fprintf (stderr, _("%s: unrecognized option `%c%s'\n"), argv[0], argv[optind][0], nextchar); } nextchar = (char *) ""; optind++; optopt = 0; return '?'; } } /* Look at and handle the next short option-character. */ { char c = *nextchar++; char *temp = my_index (optstring, c); /* Increment `optind' when we start to process its last character. */ if (*nextchar == '\0') ++optind; if (temp == NULL || c == ':') { if (opterr) { if (posixly_correct) /* 1003.2 specifies the format of this message. */ fprintf (stderr, _("%s: illegal option -- %c\n"), argv[0], c); else fprintf (stderr, _("%s: invalid option -- %c\n"), argv[0], c); } optopt = c; return '?'; } /* Convenience. Treat POSIX -W foo same as long option --foo */ if (temp[0] == 'W' && temp[1] == ';') { char *nameend; const struct option *p; const struct option *pfound = NULL; int exact = 0; int ambig = 0; int indfound = 0; int option_index; /* This is an option that requires an argument. */ if (*nextchar != '\0') { optarg = nextchar; /* If we end this ARGV-element by taking the rest as an arg, we must advance to the next element now. */ optind++; } else if (optind == argc) { if (opterr) { /* 1003.2 specifies the format of this message. */ fprintf (stderr, _("%s: option requires an argument -- %c\n"), argv[0], c); } optopt = c; if (optstring[0] == ':') c = ':'; else c = '?'; return c; } else /* We already incremented `optind' once; increment it again when taking next ARGV-elt as argument. */ optarg = argv[optind++]; /* optarg is now the argument, see if it's in the table of longopts. */ for (nextchar = nameend = optarg; *nameend && *nameend != '='; nameend++) /* Do nothing. */ ; /* Test all long options for either exact match or abbreviated matches. */ for (p = longopts, option_index = 0; p->name; p++, option_index++) if (!strncmp (p->name, nextchar, (unsigned)(nameend - nextchar))) { if ((unsigned int) (nameend - nextchar) == strlen (p->name)) { /* Exact match found. */ pfound = p; indfound = option_index; exact = 1; break; } else if (pfound == NULL) { /* First nonexact match found. */ pfound = p; indfound = option_index; } else /* Second or later nonexact match found. */ ambig = 1; } if (ambig && !exact) { if (opterr) fprintf (stderr, _("%s: option `-W %s' is ambiguous\n"), argv[0], argv[optind]); nextchar += strlen (nextchar); optind++; return '?'; } if (pfound != NULL) { option_index = indfound; if (*nameend) { /* Don't test has_arg with >, because some C compilers don't allow it to be used on enums. */ if (pfound->has_arg) optarg = nameend + 1; else { if (opterr) fprintf (stderr, _("\ %s: option `-W %s' doesn't allow an argument\n"), argv[0], pfound->name); nextchar += strlen (nextchar); return '?'; } } else if (pfound->has_arg == 1) { if (optind < argc) optarg = argv[optind++]; else { if (opterr) fprintf (stderr, _("%s: option `%s' requires an argument\n"), argv[0], argv[optind - 1]); nextchar += strlen (nextchar); return optstring[0] == ':' ? ':' : '?'; } } nextchar += strlen (nextchar); if (longind != NULL) *longind = option_index; if (pfound->flag) { *(pfound->flag) = pfound->val; return 0; } return pfound->val; } nextchar = NULL; return 'W'; /* Let the application handle it. */ } if (temp[1] == ':') { if (temp[2] == ':') { /* This is an option that accepts an argument optionally. */ if (*nextchar != '\0') { optarg = nextchar; optind++; } else optarg = NULL; nextchar = NULL; } else { /* This is an option that requires an argument. */ if (*nextchar != '\0') { optarg = nextchar; /* If we end this ARGV-element by taking the rest as an arg, we must advance to the next element now. */ optind++; } else if (optind == argc) { if (opterr) { /* 1003.2 specifies the format of this message. */ fprintf (stderr, _("%s: option requires an argument -- %c\n"), argv[0], c); } optopt = c; if (optstring[0] == ':') c = ':'; else c = '?'; } else /* We already incremented `optind' once; increment it again when taking next ARGV-elt as argument. */ optarg = argv[optind++]; nextchar = NULL; } } return c; } } int getopt ( int argc, char *const *argv, const char *optstring ) { return _getopt_internal (argc, argv, optstring, (const struct option *) 0, (int *) 0, 0); } #endif /* Not ELIDE_CODE. */