#!/bin/sh # to unpack, sh this message in an empty directory PATH=/bin:/usr/bin cat > 24048P02 <<'bigmail CUT HERE............' -iter_spcg() Sparse matrix CG method -iter_spcgne() Sparse matrix CG method for normal equations -iter_spcgs() Sparse matrix CGS method -iter_spgmres() Sparse matrix GMRES method -iter_splsqr() Sparse matrix LSQR method -iter_spmgcr() Sparse matrix MGCR method -iv_add() Add integer vectors -iv_copy() Copy integer vector -iv_dump() Dump integer vector to a stream -iv_finput() Input integer vector from a stream -iv_foutput() Output integer vector to a stream -IV_FREE() Free (deallocate) an integer vector (macro) -iv_free() Free (deallocate) integer vector (function) -iv_free_vars() Free a list of integer vectors -iv_get() Allocate and initialise an integer vector -iv_get_vars() Allocate list of integer vectors -iv_input() Input integer vector from stdin (macro) -iv_output() Output integer vector to stdout (macro) -iv_resize() Resize an integer vector -iv_resize_vars() Resize a list of integer vectors -iv_sub() Subtract integer vectors -LDLfactor() LDL^T factorisation -LDLsolve() LDL^T solver -LDLupdate() Update LDL^T factorisation -Lsolve() Solve Lx=y , L lower triangular -LTsolve() Solve L^Tx=y , L lower triangular -LUcondest() Estimate a condition number using LU factors -LUfactor() Compute LU factors with implicit scaled partial pivoting -LUsolve() Solve Ax=b using LU factors -LUTsolve() Solve A^Tx=b usng LU factors -m_add() Add matrices -makeQ() Form Q matrix for QR factorisation -makeR() Form R matrix for QR factorisation -mat2band() Extract band matrix from dense matrix -MCHfactor() Modified Cholesky factorisation - (actually factors A+D, D diagonal, instead of A) -m_copy() Copy dense matrix -m_dump() Dump matrix data structure to a stream -mem_attach_list() Adds a new family of types -mem_bytes() Notify change in memory usage (macro) -mem_bytes_list() Notify change in memory usage -mem_free_list() Frees a family of types -mem_info_bytes() Number of bytes used by a type -mem_info_numvar() Number of structures of a type -mem_info_file() Print memory info to a stream -mem_info_is_on() Is memory data being accumulated? -mem_info_on() Turns memory info system on/off -mem_is_list_attached() Is list of types attached? -mem_numvar() Notify change in number of structures allocated (macro) -mem_numvar_list() Notify change in number of structures allocated -mem_stat_dump() Prints information on registered workspace -mem_stat_free() Frees (deallocates) static workspace -mem_stat_mark() Sets mark for workspace -MEM_STAT_REG() Register static workspace (macro) -mem_stat_show_mark() Current workspace group -m_exp() Computes matrix exponential -m_finput() Input matrix from a stream -m_foutput() Output matrix to a stream -M_FREE() Free (deallocate) a matrix (macro) -m_free() Free (deallocate) matrix (function) -m_free_vars() Free a list of matrices -m_get() Allocate and initialise a matrix -m_get_vars() Allocate list of matrices -m_ident() Sets matrix to identity matrix -m_input() Input matrix from stdin (macro) -m_inverse() Invert matrix -m_load() Load matrix in MATLAB format -m_mlt() Multiplies matrices -mmtr_mlt() Computes AB^T -m_norm1() Computes ||A||_1 of a matrix -m_norm_frob() Computes the Frobenius norm of a matrix -m_norm_inf() Computes ||A||_inf of a matrix -m_ones() Set matrix to all 1's -m_output() Output matrix to stdout (macro) -m_poly() Computes a matrix polynomial -m_pow() Computes integer power of a matrix -mrand() Generates pseudo-random real number -m_rand() Randomise entries of a matrix -mrandlist() Generates array of pseudo-random numbers -m_resize() Resize matrix -m_resize_vars() Resize a list of matrices -m_save() Save matrix in MATLAB format -m_sub() Subtract matrices -m_transp() Transpose matrix -mtrm_mlt() Computes A^TB -mv_mlt() Computes Ax -mv_mltadd() Computes y <- Ax+y -m_zero() Zero a matrix -ON_ERROR() Error handler (macro) -prompter() Print prompt message to stdout -px_cols() Permute the columns of a matrix -px_copy() Copy permutation -px_dump() Dump permutation data structure to a stream -px_finput() Input permutation from a stream -px_foutput() Output permutation to a stream -PX_FREE() Free (deallocate) a permutation (macro) -px_free() Free (deallocate) permutation (function) -px_free_vars() Free a list of permutations -px_get() Allocate and initialise a permutation -px_get_vars() Allocate a list of permutations -px_ident() Sets permutation to identity -px_input() Input permutation from stdin (macro) -px_inv() Invert permutation -pxinv_vec() Computes P^Tx where P is a permutation matrix -pxinv_zvec() Computes P^Tx where P is a permutation matrix (complex) -px_mlt() Multiply permutations -px_output() Output permutation to stdout (macro) -px_resize() Resize a permutation -px_resize_vars() Resize a list of permutations -px_rows() Permute the rows of a matrix -px_sign() Returns the sign of the permutation -px_transp() Transpose a pair of entries -px_vec() Computes Px where P is a permutation matrix -px_zvec() Computes Px where P is a permutation matrix (complex) -QRCPfactor() QR factorisation with column pivoting -QRfactor() QR factorisation -QRsolve() Solve Ax=b using QR factorisation -QRTsolve() Solve A^Tx=b using QR factorisation -QRupdate() Update explicit QR factors -rot_cols() Apply Givens rotation to the columns of a matrix -rot_rows() Apply Givens rotation to the rows of a matrix -rot_vec() Apply Givens rotation to a vector -rot_zvec() Apply complex Givens rotation to a vector -schur() Compute real Schur form -schur_evals() Compute eigenvalues from the real Schur form -schur_vecs() Compute eigenvectors from the real Schur form -set_col() Set the column of a matrix to a given vector -set_err_flag() Control behaviour of ev_err() -set_row() Set the row of a matrix to a given vector -sm_mlt() Scalar-matrix multiplication -smrand() Set seed for mrand() -spBKPfactor() Sparse symmetric indefinite factorsiation -spBKPsolve() Sparse symmetric indefinite solver -spCHfactor() Sparse Cholesky factorisation -spCHsolve() Sparse Cholesky solver -spCHsymb() Symbolic sparse Cholesky factorisation - (no floating point operations) -sp_col_access() Sets up column access paths for a sparse matrix -sp_compact() Eliminates zero entries in a sparse matrix -sp_copy() Copies a sparse matrix -sp_copy2() Copies a sparse matrix into another -sp_diag_access() Sets up diagonal access paths for a sparse matrix -sp_dump() Dump sparse matrix data structure to a stream -sp_finput() Input sparse matrix from a stream -sp_foutput() Output a sparse matrix to a stream -sp_free() Free (deallocate) a sparse matrix -sp_get() Allocate and initialise a sparse matrix -sp_get_val() Get the (i,j) entry of a sparse matrix -spICHfactor() Sparse incomplete Cholesky factorisation -sp_input() Input a sparse matrix form stdin -spLUfactor() Sparse LU factorisation using partial pivoting -spLUsolve() Solves Ax=b using sparse LU factors -spLUTsolve() Solves A^Tx=b using sparse LU factors -sp_mv_mlt() Computes Ax for sparse A -sp_output() Outputs a sparse matrix to a stream (macro) -sp_resize() Resize a sparse matrix -sprow_add() Adds a pair of sparse rows -sprow_foutput() Output sparse row to a stream -sprow_get() Allocate and initialise a sparse row -sprow_get_idx() Get location of an entry in a sparse row -sprow_merge() Merge two sparse rows -sprow_mltadd() Sparse row vector multiply-and-add -sprow_set_val() Set an entry in a sparse row -sprow_smlt() Multiplies a sparse row by a scalar -sprow_sub() Subtracts a sparse row from another -sprow_xpd() Expand a sparse row -sp_set_val() Set the (i,j) entry of a sparse matrix -sp_vm_mlt() Compute x^TA for sparse A -sp_zero() Zero (but do not remove) all entries of a sparse matrix -svd() Compute the SVD of a matrix -sv_mlt() Scalar-vector multiply -symmeig() Compute eigenvalues/vectors of a symmetric matrix -tracecatch() Catch and re-raise errors (macro) -trieig() Compute eigenvalues/vectors of a symmetric tridiagonal matrix -Usolve() Solve Ux=b where U is upper triangular -UTsolve() Solve U^Tx=b where U is upper triangular -v_add() Add vectors -v_conv() Convolution product of vectors -v_copy() Copy vector -v_dump() Dump vector data structure to a stream -v_finput() Input vector from a stream -v_foutput() Output vector to a stream -V_FREE() Free (deallocate) a vector (macro) -v_free() Free (deallocate) vector (function) -v_free_vars() Free a list of vectors -v_get() Allocate and initialise a vector -v_get_vars() Allocate list of vectors -v_input() Input vector from stdin (macro) -v_lincomb() Compute sum of a_i x_i for an array of vectors -v_linlist() Compute sum of a_i x_i for a list of vectors -v_map() Apply function componentwise to a vector -v_max() Computes max vector entry and index -v_min() Computes min vector entry and index -v_mltadd() Computes y <- alpha*x+y for vectors x , y -vm_mlt() Computes x^TA -vm_mltadd() Computes y^T <- y^T+x^TA -v_norm1() Computes ||x||_1 for a vector -v_norm2() Computes ||x||_2 (the Euclidean norm) of a vector -v_norm_inf() Computes ||x||_inf for a vector -v_ones() Set vector to all 1's -v_output() Output vector to stdout (macro) -v_pconv() Periodic convolution of two vectors -v_rand() Randomise entries of a vector -v_resize() Resize a vector -v_resize_vars() Resize a list of vectors -v_save() Save a vector in MATLAB format -v_slash() Computes componentwise ratio of vectors -v_sort() Sorts vector components -v_star() Componentwise vector product -v_sub() Subtract two vectors -v_sum() Sum of components of a vector -v_zero() Zero a vector -zabs() Complex absolute value (modulus) -zadd() Add complex numbers -zconj() Conjugate complex number -zdiv() Divide complex numbers -zexp() Complex exponential -z_finput() Read complex number from file or stream -z_foutput() Prints complex number to file or stream -zgivens() Compute complex Givens' rotation -zhhtrcols() Apply Householder transformation: PA (complex) -zhhtrrows() Apply Householder transformation: AP (complex) -zhhtrvec() Apply Householder transformation: Px (complex) -zhhvec() Compute Householder transformation -zin_prod() Complex inner product -z_input() Read complex number from stdin -zinv() Computes 1/z (complex) -zLAsolve() Solve L^*x=b , L complex lower triangular -zlog() Complex logarithm -zLsolve() Solve Lx=b , L complex lower triangular -zLUAsolve() Solve A^*x=b using complex LU factorisation - (A^* - adjoint of A, A is complex) -zLUcondest() Complex LU condition estimate -zLUfactor() Complex LU factorisation -zLUsolve() Solve Ax=b using complex LU factorisation -zm_add() Add complex matrices -zm_adjoint() Computes adjoint of complex matrix -zmake() Construct complex number from real and imaginary parts -zmakeQ() Construct Q matrix for complex QR -zmakeR() Construct R matrix for complex QR -zmam_mlt() Computes A^*B (complex) -zm_dump() Dump complex matrix to stream -zm_finput() Input complex matrix from stream -ZM_FREE() Free (deallocate) complex matrix (macro) -zm_free() Free (deallocate) complex matrix (function) -zm_free_vars() Free a list of complex matrices -zm_get() Allocate complex matrix -zm_get_vars() Allocate a list of complex matrices -zm_input() Input complex matrix from stdin -zm_inverse() Compute inverse of complex matrix -zm_load() Load complex matrix in MATLAB format -zmlt() Multiply complex numbers -zmma_mlt() Computes AB^* (complex) -zm_mlt() Multiply complex matrices -zm_norm1() Complex matrix 1-norm -zm_norm_frob() Complex matrix Frobenius norm -zm_norm_inf() Complex matrix infinity-norm -zm_rand() Randomise complex matrix -zm_resize() Resize complex matrix -zm_resize_vars() Resize a list of complex matrices -zm_save() Save complex matrix in MATLAB format -zm_sub() Subtract complex matrices -zmv_mlt() Complex matrix-vector multiply -zmv_mltadd() Complex matrix-vector multiply and add -zm_zero() Zero complex matrix -zneg() Computes -z (complex) -z_output() Print complex number to stdout -zQRCPfactor() Complex QR factorisation with column pivoting -zQRCPsolve() Solve Ax = b using complex QR factorisation -zQRfactor() Complex QR factorisation -zQRAsolve() Solve A^*x = b using complex QR factorisation -zQRsolve() Solve Ax = b using complex QR factorisation -zrot_cols() Complex Givens' rotation of columns -zrot_rows() Complex Givens' rotation of rows -z_save() Save complex number in MATLAB format -zschur() Complex Schur factorisation -zset_col() Set column of complex matrix -zset_row() Set row of complex matrix -zsm_mlt() Complex scalar-matrix product -zsqrt() Square root z (complex) -zsub() Subtract complex numbers -zUAsolve() Solve U^*x=b , U complex upper triangular -zUsolve() Solve Ux=b , U complex upper triangular -zv_add() Add complex vectors -zv_copy() Copy complex vector -zv_dump() Dump complex vector to a stream -zv_finput() Input complex vector from a stream -ZV_FREE() Free (deallocate) complex vector (macro) -zv_free() Free (deallocate) complex vector (function) -zv_free_vars() Free a list of complex vectors -zv_get() Allocate complex vector -zv_get_vars() Allocate a list of complex vectors -zv_input() Input complex vector from a stdin -zv_lincomb() Compute sum of a_i x_i for an array of vectors -zv_linlist() Compute sum of a_i x_i for a list of vectors -zv_map() Apply function componentwise to a complex vector -zv_mlt() Complex scalar-vector product -zv_mltadd() Complex scalar-vector multiply and add -zvm_mlt() Computes A^*x (complex) -zvm_mltadd() Computes A^*x+y (complex) -zv_norm1() Complex vector 1-norm vnorm1() -zv_norm2() Complex vector 2-norm (Euclidean norm) -zv_norm_inf() Complex vector infinity- (or supremum) norm -zv_rand() Randomise complex vector -zv_resize() Resize complex vector -zv_resize_vars() Resize a list of complex vectors -zv_save() Save complex vector in MATLAB format -zv_slash() Componentwise ratio of complex vectors -zv_star() Componentwise product of complex vectors -zv_sub() Subtract complex vectors -zv_sum() Sum of components of a complex vector -zv_zero() Zero complex vector - - - - Low level routines - - - Function Description - -__add__() Add arrays -__ip__() Inner product of arrays -MEM_COPY() Copy memory (macro) -MEM_ZERO() Zero memory (macro) -__mltadd__() Forms x+ alpha*y for arrays -__smlt__() Scalar-vector multiplication for arrays -__sub__() Subtract an array from another -__zadd__() Add complex arrays -__zconj__() Conjugate complex array -__zero__() Zero an array -__zip__() Complex inner product of arrays -__zmlt__() Complex array scalar product -__zmltadd__() Complex array saxpy -__zsub__() Subtract complex arrays -__zzero__() Zero a complex array - - //GO.SYSIN DD DOC/fnindex.txt echo DOC/tutorial.txt 1>&2 sed >DOC/tutorial.txt <<'//GO.SYSIN DD DOC/tutorial.txt' 's/^-//' - - - MESCHACH VERSION 1.2A - --------------------- - - - TUTORIAL - ======== - - - In this manual the basic data structures are introduced, and some of the -more basic operations are illustrated. Then some examples of how to use -the data structures and procedures to solve some simple problems are given. -The first example program is a simple 4th order Runge-Kutta solver for -ordinary differential equations. The second is a general least squares -equation solver for over-determined equations. The third example -illustrates how to solve a problem involving sparse matrices. These -examples illustrate the use of matrices, matrix factorisations and solving -systems of linear equations. The examples described in this manual are -implemented in tutorial.c. - - While the description of each aspect of the system is brief and far from -comprehensive, the aim is to show the different aspects of how to set up -programs and routines and how these work in practice, which includes I/O -and error-handling issues. - - - -1. THE DATA STRUCTURES AND SOME BASIC OPERATIONS - - The three main data structures are those describing vectors, matrices -and permutations. These have been used to create data structures for -simplex tableaus for linear programming, and used with data structures for -sparse matrices etc. To use the system reliably, you should always use -pointers to these data structures and use library routines to do all the -necessary initialisation. - - In fact, for the operations that involve memory management (creation, -destruction and resizing), it is essential that you use the routines -provided. - - For example, to create a matrix A of size 34 , a vector x of dimension -10, and a permutation p of size 10, use the following code: - - - #include "matrix.h" - .............. - main() - { - MAT *A; - VEC *x; - PERM *p; - .......... - A = m_get(3,4); - x = v_get(10); - p = px_get(10); - .......... - } - - - This initialises these data structures to have the given size. The -matrix A and the vector x are initially all zero, while p is initially the -identity permutation. - - They can be disposed of by calling M_FREE(A), V_FREE(x) and PX_FREE(p) -respectively if you need to re-use the memory for something else. The -elements of each data structure can be accessed directly using the members -(or fields) of the corresponding structures. For example the (i,j) -component of A is accessed by A->me[i][j], x_i by x->ve[i] and p_i by -p->pe[i]. - - Their sizes are also directly accessible: A->m and A->n are the number -of rows and columns of A respectively, x->dim is the dimension of x , and -size of p is p->size. - - Note that the indexes are zero relative just as they are in ordinary C, -so that the index i in x->ve[i] can range from 0 to x->dim -1 . Thus the -total number of entries of a vector is exactly x->dim. - - While this alone is sufficient to allow a programmer to do any desired -operation with vectors and matrices it is neither convenient for the -programmer, nor efficient use of the CPU. A whole library has been -implemented to reduce the burden on the programmer in implementing -algorithms with vectors and matrices. For instance, to copy a vector from -x to y it is sufficient to write y = v_copy(x,VNULL). The VNULL is the -NULL vector, and usually tells the routine called to create a vector for -output. - - Thus, the v_copy function will create a vector which has the same size -as x and all the components are equal to those of x. If y has already -been created then you can write y = v_copy(x,y); in general, writing -``v_copy(x,y);'' is not enough! If y is NULL, then it is created (to have -the correct size, i.e. the same size as x), and if it is the wrong size, -then it is resized to have the correct size (i.e. same size as x). Note -that for all the following functions, the output value is returned, even if -you have a non-NULL value as the output argument. This is the standard -across the entire library. - - Addition, subtraction and scalar multiples of vectors can be computed by -calls to library routines: v_add(x,y,out), v_sub(x,y,out), sv_mlt(s,x,out) -where x and y are input vectors (with data type VEC *), out is the output -vector (same data type) and s is a double precision number (data type -double). There is also a special combination routine, which computes -out=v_1+s,v_2 in a single routine: v_mltadd(v1,v2,s,out). This is not only -extremely useful, it is also more efficient than using the scalar-vector -multiply and vector addition routines separately. - - Inner products can be computed directly: in_prod(x,y) returns the inner -product of x and y. Note that extended precision evaluation is not -guaranteed. The standard installation options uses double precision -operations throughout the library. - - Equivalent operations can be performed on matrices: m_add(A,B,C) which -returns C=A+B , and sm_mlt(s,A,C) which returns C=sA . The data types of -A, B and C are all MAT *, while that of s is type double as before. The -matrix NULL is called MNULL. - - Multiplying matrices and vectors can be done by a single function call: -mv_mlt(A,x,out) returns out=A*x while vm_mlt(A,x,out) returns out=A^T*x , or -equivalently, out^T=x^T*A . Note that there is no distinction between row -and column vectors unlike certain interactive environments such as MATLAB -or MATCALC. - - Permutations are also an essential part of the package. Vectors can be -permuted by using px_vec(p,x,p_x), rows and columns of matrices can be -permuted by using px_rows(p,A,p_A), px_cols(p,A,A_p), and permutations can -be multiplied using px_mlt(p1,p2,p1_p2) and inverted using px_inv(p,p_inv). -The NULL permutation is called PXNULL. - - There are also utility routines to initialise or re-initialise these -data structures: v_zero(x), m_zero(A), m_ident(A) (which sets A=I of the -correct size), v_rand(x), m_rand(A) which sets the entries of x and A -respectively to be randomly and uniformly selected between zero and one, -and px_ident(p) which sets p to be an identity permutation. - - Input and output are accomplished by library routines v_input(x), -m_input(A), and px_input(p). If a null object is passed to any of these -input routines, all data will be obtained from the input file, which is -stdin. If input is taken from a keyboard then the user will be prompted -for all the data items needed; if input is taken from a file, then the -input will have to be of the same format as that produced by the output -routines, which are: v_output(x), m_output(A) and px_output(p). This -output is both human and machine readable! - - If you wish to send the data to a file other than the standard output -device stdout, or receive input from a file or device other than the -standard input device stdin, take the appropriate routine above, use the -``foutpout'' suffix instead of just ``output'', and add a file pointer as -the first argument. For example, to send a matrix A to a file called -``fred'', use the following: - - - #include "matrix.h" - ............. - main() - { - FILE *fp; - MAT *A; - ............. - fp = fopen("fred","w"); - m_foutput(fp,A); - ............. - } - - - These input routines allow for the presence of comments in the data. A -comment in the input starts with a ``hash'' character ``#'', and continues -to the end of the line. For example, the following is valid input for a -3-dimensional vector: - - # The initial vector must not be zero - # x = - Vector: dim: 3 - -7 0 3 - - - For general input/output which conforms to this format, allowing -comments in the input files, use the input() and finput() macros. These -are used to print out a prompt message if stdin is a terminal (or ``tty'' -in Unix jargon), and to skip over any comments if input is from a -non-interactive device. An example of the usage of these macros is: - - input("Input number of steps: ","%d",&steps); - fp = stdin; - finput(fp,"Input number of steps: ","%d",&steps); - fp = fopen("fred","r"); - finput(fp,"Input number of steps: ","%d",&steps); - -The "%d" is one of the format specifiers which are used in fscanf(); the -last argument is the pointer to the variable (unless the variable is a -string) just as for scanf() and fscanf(). The first two macro calls read -input from stdin, the last from the file fred. If, in the first two calls, -stdin is a keyboard (a ``tty'' in Unix jargon) then the prompt string - "Input number of steps: " -is printed out on the terminal. - - - The second part of the library contains routines for various -factorisation methods. To use it put - - #include "matrix2.h" - -at the beginning of your program. It contains factorisation and solution -routines for LU, Cholesky and QR-factorisation methods, as well as update -routines for Cholesky and QR factorisations. Supporting these are a number -of Householder transformation and Givens' rotation routines. Also there is -a routine for generating the Q matrix for a QR-factorisation, if it is -needed explicitly, as it often is. -There are routines for band factorisation and solution for LU and LDL^T -factorisations. - -For using complex numbers, vectors and matrices include - - #include "zmatrix.h" - -for using the basic routines, and - - #include "zmatrix2.h" - -for the complex matrix factorisation routines. The zmatrix2.h file -includes matrix.h and zmatrix.h so you don't need these files included -together. - -For using the sparse matrix routines in the library you need to put - - #include "sparse.h" - -or, if you use any sparse factorisation routines, - - #include "sparse2.h" - -at the beginning of your file. The routines contained in the library -include routines for creating, destroying, initialising and updating sparse -matrices, and also routines for sparse matrix-dense vector multiplication, -sparse LU factorisation and sparse Cholesky factorisation. - -For using the iterative routines you need to use - - #include "iter.h" - -This includes the sparse.h and matrix.h file. -There are also routines for applying iterative methods such as -pre-conditioned conjugate gradient methods to sparse matrices. - - And if you use the standard maths library (sin(), cos(), tan(), exp(), -log(), sqrt(), acos() etc.) don't forget to include the standard -mathematics header: - - #include - -This file is not automatically included by any of the Meschach -header files. - - - -2. HOW TO MANAGE MEMORY - - Unlike many other numerical libraries, Meschach allows you to allocate, -deallocate and resize the vectors, matrices and permutations that you are -using. To gain maximum benefit from this it is sometimes necessary to -think a little about where memory is allocated and deallocated. There are -two reasons for this. - - Memory allocation, deallocation and resizing takes a significant amount -of time compared with (say) vector operations, so it should not be done too -frequently. Allocating memory but not deallocating it means that it cannot -be used by any other data structure. Data structures that are no longer -needed should be explicitly deallocated, or kept as static variables for -later use. Unlike other interpreted systems (such as Lisp) there is no -implicit ``garbage collection'' of no-longer-used memory. - - There are three main strategies that are recommended for deciding how to -allocate, deallocate and resize objects. These are ``no deallocation'' -which is really only useful for demonstration programs, ``allocate and -deallocate'' which minimises overall memory requirements at the expense of -speed, and ``resize on demand'' which is useful for routines that are -called repeatedly. A new technique for static workspace arrays is to -``register workspace variables''. - - -2.1 NO DEALLOCATION - - This is the strategy of allocating but never deallocating data -structures. This is only useful for demonstration programs run with small -to medium size data structures. For example, there could be a line - - QR = m_copy(A,MNULL); /* allocate memory for QR */ - -to allocate the memory, but without the call M_FREE(QR); in it. This can -be acceptable if QR = m_copy(A,MNULL) is only executed once, and so the -allocated memory never needs to be explicitly deallocated. - - This would not be acceptable if QR = m_copy(A,MNULL) occurred inside a -for loop. If this were so, then memory would be ``lost'' as far as the -program is concerned until there was insufficient space for allocating the -next matrix for QR. The next subsection shows how to avoid this. - - -2.2 ALLOCATE AND DEALLOCATE - - This is the most straightforward way of ensuring that memory is not -lost. With the example of allocating QR it would work like this: - - for ( ... ; ... ; ... ) - { - QR = m_copy(A,MNULL); /* allocate memory for QR */ - /* could have been allocated by m_get() */ - /* use QR */ - ...... - ...... - /* no longer need QR for this cycle */ - M_FREE(QR); /* deallocate QR so memory can be reused */ - } - - The allocate and deallocate statements could also have come at the -beginning and end of a function or procedure, so that when the function -returns, all the memory that the function has allocated has been -deallocated. - - This is most suitable for functions or sections of code that are called -repeatedly but involve fairly extensive calculations (at least a -matrix-matrix multiply, or solving a system of equations). - - -2.3 RESIZE ON DEMAND - - This technique reduces the time involved in memory allocation for code -that is repeatedly called or used, especially where the same size matrix or -vector is needed. For example, the vectors v1, v2, etc. in the -Runge-Kutta routine rk4() are allocated according to this strategy: - - rk4(...,x,...) - { - static VEC *v1=VNULL, *v2=VNULL, *v3=VNULL, *v4=VNULL, *temp=VNULL; - ....... - v1 = v_resize(v1,x->dim); - v2 = v_resize(v2,x->dim); - v3 = v_resize(v3,x->dim); - v4 = v_resize(v4,x->dim); - temp = v_resize(temp,x->dim); - ....... - } - - The intention is that the rk4() routine is called repeatedly with the -same size x vector. It then doesn't make as much sense to allocate v1, v2 -etc. whenever the function is called. Instead, v_resize() only performs -memory allocation if the memory already allocated to v1, v2 etc. is smaller -than x->dim. - - The vectors v1, v2 etc. are declared to be static to ensure that their -values are not lost between function calls. Variables that are declared -static are set to NULL or zero by default. So the declaration of v1, v2, -etc., could be - - static VEC *v1, *v2, *v3, *v4, *temp; - - This strategy of resizing static workspace variables is not so useful if -the object being allocated is extremely large. The previous ``allocate and -deallocate'' strategy is much more efficient for memory in those -circumstances. However, the following section shows how to get the best of -both worlds. - - -2.4 REGISTRATION OF WORKSPACE - - From version 1.2 onwards, workspace variables can be registered so that -the memory they reference can be freed up on demand. To do this, the -function containing the static workspace variables has to include calls to -MEM_STAT_REG(var,type) where var is a pointer to a Meschach data type (such -as VEC or MAT). This call should be placed after the call to the -appropriate resize function. The type parameter should be a TYPE_... macro -where the ``...'' is the name of a Meschach type such as VEC or MAT. For -example, - - rk4(...,x,...) - { - static VEC *v1, *v2, *v3, *v4, *temp; - ....... - v1 = v_resize(v1,x->dim); - MEM_STAT_REG(v1,TYPE_VEC); - v2 = v_resize(v2,x->dim); - MEM_STAT_REG(v2,TYPE_VEC); - ...... - } - -Normally, these registered workspace variables remain allocated. However, -to implement the ``deallocate on exit'' approach, use the following code: - - ...... - mem_stat_mark(1); - rk4(...,x,...) - mem_stat_free(1); - ...... - - To keep the workspace vectors allocated for the duration of a loop, but -then deallocated, use - - ...... - mem_stat_mark(1); - for (i = 0; i < N; i++ ) - rk4(...,x,...); - mem_stat_free(1); - ...... - -The number used in the mem_stat_mark() and mem_stat_free() calls is the -workspace group number. The call mem_stat_mark(1) designates 1 as the -current workspace group number; the call mem_stat_free(1) deallocates (and -sets to NULL) all static workspace variables registered as belonging to -workspace group 1. - - - -3. SIMPLE VECTOR OPERATIONS: AN RK4 ROUTINE - - The main purpose of this example is to show how to deal with vectors and -to compute linear combinations. - - The problem here is to implement the standard 4th order Runge-Kutta -method for the ODE - - x'=f(t,x), x(t_0)=x_0 - -for x(t_i), i=1,2,3, where t_i=t_0+i*h and h is the step size. - - The formulae for the 4th order Runge-Kutta method are: - - x_i+1 = x_i+ h/6*(v_1+2*v_2+2*v_3+v_4), -where - v_1 = f(t_i,x_i) - v_2 = f(t_i+h, x_i+h*v_1) - v_3 = f(t_i+h, x_i+h*v_2) - v_4 = f(t_i+h, x_i+h*v_3) - -where the v_i are vectors. - - The procedure for implementing this method (rk4()) will be passed (a -pointer to) the function f. The implementation of f could, in this system, -create a vector to hold the return value each time it is called. However, -such a scheme is memory intensive and the calls to the memory allocation -functions could easily dominate the time performed doing numerical -computations. So, the implementation of f will also be passed an already -allocated vector to be filled in with the appropriate values. - - The procedure rk4() will also be passed the current time t, the step -size h, and the current value for x. The time after the step will be -returned by rk4(). - -The code that does this follows. - - - #include "matrix.h" - - /* rk4 - 4th order Runge-Kutta method */ - double rk4(f,t,x,h) - double t, h; - VEC *(*f)(), *x; - { - static VEC *v1=VNULL, *v2=VNULL, *v3=VNULL, *v4=VNULL; - static VEC *temp=VNULL; - - /* do not work with NULL initial vector */ - if ( x == VNULL ) - error(E_NULL,"rk4"); - - /* ensure that v1, ..., v4, temp are of the correct size */ - v1 = v_resize(v1,x->dim); - v2 = v_resize(v2,x->dim); - v3 = v_resize(v3,x->dim); - v4 = v_resize(v4,x->dim); - temp = v_resize(temp,x->dim); - - /* register workspace variables */ - MEM_STAT_REG(v1,TYPE_VEC); - MEM_STAT_REG(v2,TYPE_VEC); - MEM_STAT_REG(v3,TYPE_VEC); - MEM_STAT_REG(v4,TYPE_VEC); - MEM_STAT_REG(temp,TYPE_VEC); - /* end of memory allocation */ - - (*f)(t,x,v1); /* most compilers allow: f(t,x,v1); */ - v_mltadd(x,v1,0.5*h,temp); /* temp = x+.5*h*v1 */ - (*f)(t+0.5*h,temp,v2); - v_mltadd(x,v2,0.5*h,temp); /* temp = x+.5*h*v2 */ - (*f)(t+0.5*h,temp,v3); - v_mltadd(x,v3,h,temp); /* temp = x+h*v3 */ - (*f)(t+h,temp,v4); - - /* now add: v1+2*v2+2*v3+v4 */ - v_copy(v1,temp); /* temp = v1 */ - v_mltadd(temp,v2,2.0,temp); /* temp = v1+2*v2 */ - v_mltadd(temp,v3,2.0,temp); /* temp = v1+2*v2+2*v3 */ - v_add(temp,v4,temp); /* temp = v1+2*v2+2*v3+v4 */ - - /* adjust x */ - v_mltadd(x,temp,h/6.0,x); /* x = x+(h/6)*temp */ - - return t+h; /* return the new time */ - } - - - Note that the last parameter of f() is where the output is placed. -Often this can be NULL in which case the appropriate data structure is -allocated and initialised. Note also that this routine can be used for -problems of arbitrary size, and the dimension of the problem is determined -directly from the data given. The vectors v_1,...,v_4 are created to have -the correct size in the lines - - .... - v1 = v_resize(v1,x->dim); - v2 = v_resize(v2,x->dim); - .... - - Here v_resize(v,dim) resizes the VEC structure v to hold a vector of -length dim. If v is initially NULL, then this creates a new vector of -dimension dim, just as v_get(dim) would do. For the above piece of code to -work correctly, v1, v2 etc., must be initialised to be NULL vectors. This -is done by the declaration - - static VEC *v1=VNULL, *v2=VNULL, *v3=VNULL, *v4=VNULL; - -or - - static VEC *v1, *v2, *v3, *v4; - -The operations of vector addition and scalar addition are really the only -vector operations that need to be performed in rk4. Vector addition is -done by v_add(v1,v2,out), where out=v1+v2, and scalar multiplication by -sv_mlt(scale,v,out), where out=scale*v. - -These can be combined into a single operation v_mltadd(v1,v2,scale,out), -where out=v1+scale*v2. As many operations in numerical mathematics involve -accumulating scalar multiples, this is an extremely useful operation, as we -can see above. For example: - - v_mltadd(x,v1,0.5*h,temp); /* temp = x+0.5*h*v1 */ - - We also need a number of ``utility'' operations. For example v_copy(in, -out) copies the vector in to out. There is also v_zero(v) to zero a vector -v. - - Here is an implementation of the function f for simple harmonic motion: - - /* f - right-hand side of ODE solver */ - VEC *f(t,x,out) - VEC *x, *out; - double t; - { - if ( x == VNULL || out == VNULL ) - error(E_NULL,"f"); - if ( x->dim != 2 || out->dim != 2 ) - error(E_SIZES,"f"); - - out->ve[0] = x->ve[1]; - out->ve[1] = - x->ve[0]; - - return out; - } - - As can be seen, most of this code is error checking code, which, of -course, makes the routine safer but a little slower. For a procedure like -f() it is probably not necessary, although then the main program would have -to perform checking to ensure that the vectors involved have the correct -size etc. The ith component of a vector x is x->ve[i], and indexing is -zero-relative (i.e., the ``first'' component is component 0). The ODE -described above is for simple harmonic motion: - - x_0'=x_1 , x_1'=-x_0 , or equivalently, x_0''+ x_0 = 0 . - - Here is the main program: - - - #include - #include "matrix.h" - - main() - { - VEC *x; - VEC *f(); - double h, t, t_fin; - double rk4(); - - input("Input initial time: ", "%lf", &t); - input("Input final time: ", "%lf", &t_fin); - x = v_get(2); /* this is the size needed by f() */ - prompter("Input initial state:\n"); x = v_input(VNULL); - input("Input step size: ", "%lf", &h); - - printf("# At time %g, the state is\n",t); - v_output(x); - while ( t < t_fin ) - { - t = rk4(f,t,x,min(h,t_fin-t)); /* new t is returned */ - printf("# At time %g, the state is\n",t); - v_output(x); - t += h; - } - } - - The initial values are entered as a vector by v_input(). If v_input() -is passed a vector, then this vector will be used to store the input, and -this vector has the size that x had on entry to v_input(). The original -values of x are also used as a prompt on input from a tty. If a NULL is -passed to v_input() then v_input() will return a vector of whatever size -the user inputs. So, to ensure that only a two-dimensional vector is used -for the initial conditions (which is what f() is expecting) we use - - x = v_get(2); x = v_input(x); - - To compile the program under Unix, if it is in a file tutorial.c: - - cc -o tutorial tutorial.c meschach.a - -or, if you have an ANSI compiler, - - cc -DANSI_C -o tutorial tutorial.c meschach.a - - Here is a sample session with the above program: - - tutorial - - Input initial time: 0 - Input final time: 1 - Input initial state: - Vector: dim: 2 - entry 0: -1 - entry 1: b - entry 0: old -1 new: 1 - entry 1: old 0 new: 0 - Input step size: 0.1 - At time 0, the state is - Vector: dim: 2 - 1 0 - At time 0.1, the state is - Vector: dim: 2 - 0.995004167 -0.0998333333 - ................. - At time 1, the state is - Vector: dim: 2 - 0.540302967 -0.841470478 - - By way of comparison, the state at t=1 for the true solution is - x_0(1)=0.5403023058 , x_1(1)=-0.8414709848 . -The ``b'' that is typed in entering the x vector allows the user to alter -previously entered components. In this case once this is done, the user is -prompted with the old values when entering the new values. The user can -also type in ``f'' for skipping over the vector's components, which are -then unchanged. If an incorrectly sized initial value vector x is given, -the error handler comes into action: - - Input initial time: 0 - Input final time: 1 - Input initial state: - Vector: dim: 3 - entry 0: 3 - entry 1: 2 - entry 2: -1 - Input step size: 0.1 - At time 0, the state is - Vector: dim: 3 - 3 2 -1 - - "tutorial.c", line 79: sizes of objects don't match in function f() - Sorry, aborting program - - The error handler prints out the error message giving the source code -file and line number as well as the function name where the error was -raised. The relevant section of f() in file tutorial.c is: - - if ( x->dim != 2 || out->dim != 2 ) - error(E_SIZES,"f"); /* line 79 */ - - - The standard routines in this system perform error checking of this -type, and also checking for undefined results such as division by zero in -the routines for solving systems of linear equations. There are also error -messages for incorrectly formatted input and end-of-file conditions. - - To round off the discussion of this program, note that we have seen -interactive input of vectors. If the input file or stream is not a tty -(e.g., a file, a pipeline or a device) then it expects the input to have -the same form as the output for each of the data structures. Each of the -input routines (v_input(), m_input(), px_input()) skips over ``comments'' -in the input data, as do the macros input() and finput(). Anything from a -`#' to the end of the line (or EOF) is considered to be a comment. For -example, the initial value problem could be set up in a file ivp.dat as: - - # Initial time - 0 - # Final time - 1 - # Solution is x(t) = (cos(t),-sin(t)) - # x(0) = - Vector: dim: 2 - 1 0 - # Step size - 0.1 - - The output of the above program with the above input (from a file) gives -essentially the same output as shown above, except that no prompts are sent -to the screen. - - - -4. USING ROUTINES FOR LISTS OF ARGUMENTS - - Some of the most common routines have variants that take a variable -number of arguments. These are the routines .._get_vars(), .._resize_vars() -and .._free_vars(). These correspond to the the basic routines .._get(), -.._resize() and .._free() respectively. Also there is the -mem_stat_reg_vars() routine which registers a list of static workspace -variables. This corresponds to mem_stat_reg_list() for a single variable. - - Here is an example of how to use these functions. This example also -uses the routine v_linlist() to compute a linear combination of vectors. -Note that the code is much more compact, but don't forget that these -``..._vars()'' routines usually need the address-of operator ``&'' and NULL -termination of the arguments to work correctly. - - - #include "matrix.h" - - /* rk4 - 4th order Runge-Kutta method */ - double rk4(f,t,x,h) - double t, h; - VEC *(*f)(), *x; - { - static VEC *v1, *v2, *v3, *v4, *temp; - - /* do not work with NULL initial vector */ - if ( x == VNULL ) - error(E_NULL,"rk4"); - - /* ensure that v1, ..., v4, temp are of the correct size */ - v_resize_vars(x->dim, &v1, &v2, &v3, &v4, &temp, NULL); - - /* register workspace variables */ - mem_stat_reg_vars(0, TYPE_VEC, &v1, &v2, &v3, &v4, &temp, NULL); - /* end of memory allocation */ - - (*f)(t,x,v1); v_mltadd(x,v1,0.5*h,temp); - (*f)(t+0.5*h,temp,v2); v_mltadd(x,v2,0.5*h,temp); - (*f)(t+0.5*h,temp,v3); v_mltadd(x,v3,h,temp); - (*f)(t+h,temp,v4); - - /* now add: temp = v1+2*v2+2*v3+v4 */ - v_linlist(temp, v1, 1.0, v2, 2.0, v3, 2.0, v4, 1.0, VNULL); - /* adjust x */ - v_mltadd(x,temp,h/6.0,x); /* x = x+(h/6)*temp */ - - return t+h; /* return the new time */ - } - - - -5. A LEAST SQUARES PROBLEM - - Here we need to use matrices and matrix factorisations (in particular, a -QR factorisation) in order to find the best linear least squares solution -to some data. Thus in order to solve the (approximate) equations - A*x = b, -where A is an m x n matrix (m > n) we really need to solve the optimisation -problem - min_x ||Ax-b||^2. - - If we write A=QR where Q is an orthogonal m x m matrix and R is an upper -triangular m x n matrix then (we use 2-norm) - - ||A*x-b||^2 = ||R*x-Q^T*b||^2 = || R_1*x - Q_1^T*b||^2 + ||Q_2^T*b||^2 - -where R_1 is an n x n upper triangular matrix. If A has full rank then R_1 -will be an invertible matrix, and the best least squares solution of A*x=b -is x= R_1^{-1}*Q_1^T*b . - - These calculations can be be done quite easily as there is a QRfactor() -function available with the system. QRfactor() is declared to have the -prototype - - MAT *QRfactor(MAT *A, VEC *diag); - - The matrix A is overwritten with the factorisation of A ``in compact -form''; that is, while the upper triangular part of A is indeed the R -matrix described above, the Q matrix is stored as a collection of -Householder vectors in the strictly lower triangular part of A and in the -diag vector. The QRsolve() function knows and uses this compact form and -solves Q*R*x=b with the call QRsolve(A,diag,b,x), which also returns x. - - Here is the code to obtain the matrix A, perform the QR factorisation, -obtain the data vector b, solve for x, and determine what the norm of the -errors ( ||Ax-b||_2 ) is. - - - #include "matrix2.h" - - main() - { - MAT *A, *QR; - VEC *b, *x, *diag; - - /* read in A matrix */ - printf("Input A matrix:"); - - A = m_input(MNULL); /* A has whatever size is input */ - - if ( A->m < A->n ) - { - printf("Need m >= n to obtain least squares fit"); - exit(0); - } - printf("# A ="); m_output(A); - diag = v_get(A->m); - - /* QR is to be the QR factorisation of A */ - QR = m_copy(A,MNULL); - QRfactor(QR,diag); - - /* read in b vector */ - printf("Input b vector:"); - b = v_get(A->m); - b = v_input(b); - printf("# b ="); v_output(b); - - /* solve for x */ - x = QRsolve(QR,diag,b,VNULL); - printf("Vector of best fit parameters is"); - v_output(x); - - /* ... and work out norm of errors... */ - printf("||A*x-b|| = %g\n", - v_norm2(v_sub(mv_mlt(A,x,VNULL),b,VNULL))); - } - - Note that as well as the usual memory allocation functions like m_get(), -the I/O functions like m_input() and m_output(), and the -factorise-and-solve functions QRfactor() and QRsolve(), there are also -functions for matrix-vector multiplication: - mv_mlt(MAT *A, VEC *x, VEC *out) -and also vector-matrix multiplication (with the vector on the left): - vm_mlt(MAT *A, VEC *x, VEC *out), -with out=x^T A. There are also functions to perform matrix arithmetic - -matrix addition m_add(), matrix-scalar multiplication sm_mlt(), -matrix-matrix multiplication m_mlt(). - - Several different sorts of matrix factorisation are supported: LU -factorisation (also known as Gaussian elimination) with partial pivoting, -by LUfactor() and LUsolve(). Other factorisation methods include Cholesky -factorisation CHfactor() and CHsolve(), and QR factorisation with column -pivoting QRCPfactor(). - - Pivoting involve permutations which have their own PERM data structure. -Permutations can be created by px_get(), read and written by px_input() and -px_output(), multiplied by px_mlt(), inverted by px_inv() and applied to -vectors by px_vec(). - -The above program can be put into a file leastsq.c and compiled under Unix -using - - cc -o leastsq leastsq.c meschach.a -lm - -A sample session using leastsq follows: - - - Input A matrix: - Matrix: rows cols:5 3 - row 0: - entry (0,0): 3 - entry (0,1): -1 - entry (0,2): 2 - Continue: - row 1: - entry (1,0): 2 - entry (1,1): -1 - entry (1,2): 1 - Continue: n - row 1: - entry (1,0): old 2 new: 2 - entry (1,1): old -1 new: -1 - entry (1,2): old 1 new: 1.2 - Continue: - row 2: - entry (2,0): old 0 new: 2.5 - .... - .... (Data entry) - .... - # A = - Matrix: 5 by 3 - row 0: 3 -1 2 - row 1: 2 -1 1.2 - row 2: 2.5 1 -1.5 - row 3: 3 1 1 - row 4: -1 1 -2.2 - Input b vector: - entry 0: old 0 new: 5 - entry 1: old 0 new: 3 - entry 2: old 0 new: 2 - entry 3: old 0 new: 4 - entry 4: old 0 new: 6 - # b = - Vector: dim: 5 - 5 3 2 4 6 - Vector of best fit parameters is - Vector: dim: 3 - 1.47241555 -0.402817858 -1.14411815 - ||A*x-b|| = 6.78938 - - - The Q matrix can be obtained explicitly by the routine makeQ(). The Q -matrix can then be used to obtain an orthogonal basis for the range of A . -An orthogonal basis for the null space of A can be obtained by finding the -QR-factorisation of A^T . - - - -6. A SPARSE MATRIX EXAMPLE - - To illustrate the sparse matrix routines, consider the problem of -solving Poisson's equation on a square using finite differences, and -incomplete Cholesky factorisation. The actual equations to solve are - - u_{i,j+1} + u_{i,j-1} + u_{i+1,j} + u_{i-1,j} - 4*u_{i,j} = - h^2*f(x_i,y_j), for i,j=1,...,N - -where u_{0,j} = u_{i,0} = u_{N+1,j} = u_{i,N+1} = 0 for i,j=1,...,N and h -is the common distance between grid points. - - The first task is to set up the matrix describing this system of linear -equations. The next is to set up the right-hand side. The third is to -form the incomplete Cholesky factorisation of this matrix, and finally to -use the sparse matrix conjugate gradient routine with the incomplete -Cholesky factorisation as preconditioner. - - Setting up the matrix and right-hand side can be done by the following -code: - - - #define N 100 - #define index(i,j) (N*((i)-1)+(j)-1) - ...... - A = sp_get(N*N,N*N,5); - b = v_get(N*N); - h = 1.0/(N+1); /* for a unit square */ - ...... - - for ( i = 1; i <= N; i++ ) - for ( j = 1; j <= N; j++ ) - { - if ( i < N ) - sp_set_val(A,index(i,j),index(i+1,j),-1.0); - if ( i > 1 ) - sp_set_val(A,index(i,j),index(i-1,j),-1.0); - if ( j < N ) - sp_set_val(A,index(i,j),index(i,j+1),-1.0); - if ( j > 1 ) - sp_set_val(A,index(i,j),index(i,j-1),-1.0); - sp_set_val(A,index(i,j),index(i,j),4.0); - b->ve[index(i,j)] = -h*h*f(h*i,h*j); - } - - Once the matrix and right-hand side are set up, the next task is to -compute the sparse incomplete Cholesky factorisation of A. This must be -done in a different matrix, so A must be copied. - - LLT = sp_copy(A); - spICHfactor(LLT); - -Now when that is done, the remainder is easy: - - out = v_get(A->m); - ...... - iter_spcg(A,LLT,b,1e-6,out,1000,&num_steps); - printf("Number of iterations = %d\n",num_steps); - ...... - -and the output can be used in whatever way desired. - - For graphical output of the results, the solution vector can be copied -into a square matrix, which is then saved in MATLAB format using m_save(), -and graphical output can be produced by MATLAB. - - - -7. HOW DO I ....? - - For the convenience of the user, here a number of common tasks that -people need to perform frequently, and how to perform the computations -using Meschach. - - -7.1 .... SOLVE A SYSTEM OF LINEAR EQUATIONS ? - - If you wish to solve Ax=b for x given A and b (without destroying A), -then the following code will do this: - - VEC *x, *b; - MAT *A, *LU; - PERM *pivot; - ...... - LU = m_get(A->m,A->n); - LU = m_copy(A,LU); - pivot = px_get(A->m); - LUfactor(LU,pivot); - /* set values of b here */ - x = LUsolve(LU,pivot,b,VNULL); - - -7.2 .... SOLVE A LEAST-SQUARES PROBLEM ? - - To minimise ||Ax-b||_2^2 = sum_i ((Ax)_i-b_i)^2, the most reliable -method is based on the QR-factorisation. The following code performs this -calculation assuming that A is m x n with m > n : - - MAT *A, *QR; - VEC *diag, *b, *x; - ...... - QR = m_get(A->m,A->n); - QR = m_copy(A,QR); - diag = v_get(A->n); - QRfactor(QR,diag); - /* set values of b here */ - x = QRsolve(QR,diag,b,x); - - -7.3 .... FIND ALL THE EIGENVALUES (AND EIGENVECTORS) OF A GENERAL MATRIX ? - - The best method is based on the Schur decomposition. For symmetric -matrices, the eigenvalues and eigenvectors can be computed by a single call -to symmeig(). For non-symmetric matrices, the situation is more complex -and the problem of finding eigenvalues and eigenvectors can become quite -ill-conditioned. Provided the problem is not too ill-conditioned, the -following code should give accurate results: - - - /* A is the matrix whose eigenvalues and eigenvectors are sought */ - MAT *A, *T, *Q, *X_re, *X_im; - VEC *evals_re, *evals_im; - ...... - Q = m_get(A->m,A->n); - T = m_copy(A,MNULL); - - /* compute Schur form: A = Q*T*Q^T */ - schur(T,Q); - - /* extract eigenvalues */ - evals_re = v_get(A->m); - evals_im = v_get(A->m); - schur_evals(T,evals_re,evals_im); - - /* Q not needed for eiegenvalues */ - X_re = m_get(A->m,A->n); - X_im = m_get(A->m,A->n); - schur_vecs(T,Q,X_re,X_im); - /* k'th eigenvector is k'th column of (X_re + i*X_im) */ - - - -7.4 .... SOLVE A LARGE, SPARSE, POSITIVE DEFINITE SYSTEM OF EQUATIONS ? - - An example of a large, sparse, positive definite matrix is the matrix -obtained from a finite-difference approximation of the Laplacian operator. -If an explicit representation of such a matrix is available, then the -following code is suggested as a reasonable way of computing solutions: - - - /* A*x == b is the system to be solved */ - SPMAT *A, *LLT; - VEC *x, *b; - int num_steps; - ...... - /* set up A and b */ - ...... - x = m_get(A->m); - LLT = sp_copy(A); - - /* preconditioning using the incomplete Cholesky factorisation */ - spICHfactor(LLT); - - /* now use pre-conditioned conjugate gradients */ - x = iter_spcg(A,LLT,b,1e-7,x,1000,&num_steps); - /* solution computed to give a relative residual of 10^-7 */ - - - If explicitly storing such a matrix takes up too much memory, then if -you can write a routine to perform the calculation of A*x for any given x , -the following code may be more suitable (if slower): - - - VEC *mult_routine(user_def,x,out) - void *user_def; - VEC *x, *out; - { - /* compute out = A*x */ - ...... - return out; - } - - - main() - { - ITER *ip; - VEC *x, *b; - ...... - b = v_get(BIG_DIM); /* right-hand side */ - x = v_get(BIG_DIM); /* solution */ - - /* set up b */ - ...... - ip = iter_get(b->dim, x->dim); - ip->b = v_copy(b,ip->b); - ip->info = NULL; /* if you don't want information - about solution process */ - v_zero(ip->x); /* initial guess is zero */ - iter_Ax(ip,mult_routine,user_def); - iter_cg(ip); - printf("# Solution is:\n"); v_output(ip->x); - ...... - ITER_FREE(ip); /* destroy ip */ - } - - The user_def argument is for a pointer to a user-defined structure -(possibly NULL, if you don't need this) so that you can write a common -function for handling a large number of different circumstances. - - - -8. MORE ADVANCED TOPICS - - Read this if you are interested in using Meschach library as a base for -applications. As an example we show how to implement a new type for 3 -dimensional matrices and incorporate this new type into the Meschach -system. Usually this part of Meschach is transparent to a user. But a more -advanced user can take advantage of these routines. We do not describe -the routines in detail here, but we want to give a rather broad picture of -what can be done. By the system we mainly mean the system of delivering -information on the number of bytes of allocated memory and routines for -deallocating static variables by mem_stat_... routines. - - First we introduce a concept of a list of types. By a list of types we -mean a set of different types with corresponding routines for creating -these types, destroying and resizing them. Each type list has a number. -The list 0 is a list of standard Meschach types such as MAT or VEC. Other -lists can be defined by a user or a application (based on Meschach). The -user can attach his/her own list to the system by the routine -mem_attach_list(). Sometimes it is worth checking if a list number is -already used by another application. It can be done by -mem_is_list_attached(ls_num), which returns TRUE if the number ls_num -is used. And such a list can be removed from the system by -mem_free_list(ls_num) if necessary. - - We describe arguments required by mem_attach_list(). The prototype of -this function is as follow - - int mem_attach_list(int ls_num, int ntypes, char *type_names[], - int (*free_funcs[])(), MEM_ARRAY sum[]); - -where the structure MEM_ARRAY has two members: "bytes" of type long and -"numvar" of type int. The frst argument is the list number. Note that you -cannot overwrite another list. To do this remove first the old list (by -mem_free_list()) or choose another number. The next argument is the number -of types which are on the list. This number cannot be changed during -running a program. The third argument is an array containing the names of -types (these are character strings). The fourth one is an array of -functions deallocating variables of the corresponding type. And the last -argument is the local array where information about the number of bytes of -allocated/deallocated memory (member bytes) and the number of allocated -variables (member numvar) are gathered. The functions which send -information to this array are mem_bytes_list() and mem_numvar_list(). - - -Example: The routines described here are in the file tutadv.c. -Firstly we define some macros and a type for 3 dimensional matrices. - -#include "matrix.h" -#define M3D_LIST 3 /* list number */ -#define TYPE_MAT3D 0 /* the number of a type */ -/* type for 3 dimensional matrices */ -typedef struct { - int l,m,n; /* actual dimensions */ - int max_l, max_m, max_n; /* maximal dimensions */ - Real ***me; /* pointer to matrix elements */ - /* we do not consider segmented memory */ - Real *base, **me2d; /* me and me2d are additional pointers - to base */ -} MAT3D; - - -Now we need two routines: one for allocating memory for 3 dimensional -matrices and the other for deallocating it. It can be useful to have a -routine for resizing 3 dimensional matrices but we do not use it here. -Note the use of mem_bytes_list() and mem_numvar_list() to notify the change -in the number of structures and bytes in use. - -/* function for creating a variable of MAT3D type */ - -MAT3D *m3d_get(l,m,n) -int l,m,n; -{ - MAT3D *mat; - .... - /* alocate memory for structure */ - if ((mat = NEW(MAT3D)) == (MAT3D *)NULL) - error(E_MEM,"m3d_get"); - else if (mem_info_is_on()) { - /* record how many bytes are allocated to structure */ - mem_bytes_list(TYPE_MAT3D,0,sizeof(MAT3D),M3D_LIST); - /* record a new allocated variable */ - mem_numvar_list(TYPE_MAT3D,1,M3D_LIST); - } - .... - /* allocate memory for 3D array */ - if ((mat->base = NEW_A(l*m*n,Real)) == (Real *)NULL) - error(E_MEM,"m3d_get"); - else if (mem_info_is_on()) - mem_bytes_list(TYPE_MAT3D,0,l*m*n*sizeof(Real),M3D_LIST); - .... - return mat; -} - -/* deallocate a variable of type MAT3D */ - -int m3d_free(mat) -MAT3D *mat; -{ - /* do not try to deallocate the NULL pointer */ - if (mat == (MAT3D *)NULL) - return -1; - .... - /* first deallocate base */ - if (mat->base != (Real *)NULL) { - if (mem_info_is_on()) - /* record how many bytes is deallocated */ - mem_bytes_list(TYPE_MAT3D,mat->max_l*mat->max_m*mat->max_n*sizeof(Real), - 0,M3D_LIST); - free((char *)mat->base); - } - .... - /* deallocate MAT3D structure */ - if (mem_info_is_on()) { - mem_bytes_list(TYPE_MAT3D,sizeof(MAT3D),0,M3D_LIST); - mem_numvar_list(TYPE_MAT3D,-1,M3D_LIST); - } - free((char *)mat); - - .... - free((char *)mat); - - return 0; -} - - -We can now create the arrays necessary for mem_attach_list(). Note that -m3d_sum can be static if it is in the same file as main(), where -mem_attach_list is called. Otherwise it must be global. - - -char *m3d_names[] = { - "MAT3D" -}; - -#define M3D_NUM (sizeof(m3d_names)/sizeof(*m3d_names)) - -int (*m3d_free_funcs[M3D_NUM])() = { - m3d_free -} - -static MEM_ARRAY m3d_sum[M3D_NUM]; - - -The last thing is to attach the list to the system. - -void main() -{ - MAT3D *M; - .... - - mem_info_on(TRUE); /* switch memory info on */ - /* attach the new list */ - mem_attach_list(M3D_LIST,M3D_NUM,m3d_names,m3d_free_funcs,m3d_sum); - .... - M = m3d_get(3,4,5); - .... - /* making use of M->me[i][j][k], where i,j,k are non-negative and - i < 3, j < 4, k < 5 */ - .... - mem_info_file(stdout,M3D_LIST); /* info on the number of allocated - bytes of memory for types - on the list M3D_LIST */ - .... - m3d_free(M); /* if M is not necessary */ - .... -} - - -We can now use the function mem_info_file() for getting information about -the number of bytes of allocated memory and number of allocated variables -of type MAT3D; mem_stat_reg_list() for registering variables of this type -and mem_stat_mark() and mem_stat_free_list() for deallocating static -variables of this type. - - - -In the similar way you can create you own list of errors and attach it to -the system. See the functions: - - int err_list_attach(int list_num, int list_len, char **err_ptr, - int warn); /* for attaching a list of errors */ - - 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 */ - -where list_num is the number of the error list, list_len is the number of -errors on the list, err_ptr is the character string explaining the error -and warn can be TRUE if this is only a warning (the program continues to -run) or it can be FALSE if it is an error (the program stops). - -The examples are the standard errors (error list 0) and warnings -(error list 1) which are in the file err.c - - - David Stewart and Zbigniew Leyk, 1993 //GO.SYSIN DD DOC/tutorial.txt mkdir MACHINES mkdir MACHINES/GCC echo MACHINES/GCC/makefile 1>&2 sed >MACHINES/GCC/makefile <<'//GO.SYSIN DD MACHINES/GCC/makefile' 's/^-//' -# -# -# Makefile for Meschach for GNU cc -# -# Copyright (C) David Stewart & Zbigniew Leyk 1993 -# -# $Id: m4,v 1.1.1.1 1999/04/14 14:16:22 borland Exp $ -# - -srcdir = . -VPATH = . - -CC = gcc - -DEFS = -DHAVE_CONFIG_H -LIBS = -lm -RANLIB = ranlib - - -CFLAGS = -O6 - - -.c.o: - $(CC) -c $(CFLAGS) $(DEFS) $< - -SHELL = /bin/sh -MES_PAK = mesch12a -TAR = tar -SHAR = stree -u -ZIP = zip -r -l - -############################### - -LIST1 = copy.o err.o matrixio.o memory.o vecop.o matop.o pxop.o \ - submat.o init.o otherio.o machine.o matlab.o ivecop.o version.o \ - meminfo.o memstat.o -LIST2 = lufactor.o bkpfacto.o chfactor.o qrfactor.o solve.o hsehldr.o \ - givens.o update.o norm.o hessen.o symmeig.o schur.o svd.o fft.o \ - mfunc.o bdfactor.o -LIST3 = sparse.o sprow.o sparseio.o spchfctr.o splufctr.o \ - spbkp.o spswap.o iter0.o itersym.o iternsym.o -ZLIST1 = zmachine.o zcopy.o zmatio.o zmemory.o zvecop.o zmatop.o znorm.o \ - zfunc.o -ZLIST2 = zlufctr.o zsolve.o zmatlab.o zhsehldr.o zqrfctr.o \ - zgivens.o zhessen.o zschur.o - -# they are no longer supported -# if you use them add oldpart to all and sparse -OLDLIST = conjgrad.o lanczos.o arnoldi.o - -ALL_LISTS = $(LIST1) $(LIST2) $(LIST3) $(ZLIST1) $(ZLIST2) $(OLDLIST) - - -HLIST = err.h iter.h machine.h matlab.h matrix.h matrix2.h \ - meminfo.h oldnames.h sparse.h sparse2.h \ - zmatrix.h zmatrix2.h - -TORTURE = torture.o sptort.o ztorture.o memtort.o itertort.o \ - mfuntort.o iotort.o - -OTHERS = dmacheps.c extras.c fmacheps.c maxint.c makefile.in \ - README configure configure.in machine.h.in copyright \ - tutorial.c tutadv.c rk4.dat ls.dat makefile - - -# Different configurations -all: part1 part2 part3 zpart1 zpart2 -basic: part1 part2 -sparse: part1 part2 part3 -complex: part1 part2 zpart1 zpart2 - - -HBASE = err.h meminfo.h machine.h matrix.h - -$(LIST1): $(HBASE) -part1: $(LIST1) - ar ru meschach.a $(LIST1); $(RANLIB) meschach.a - -$(LIST2): $(HBASE) matrix2.h -part2: $(LIST2) - ar ru meschach.a $(LIST2); $(RANLIB) - -$(LIST3): $(HBASE) sparse.h sparse2.h -part3: $(LIST3) - ar ru meschach.a $(LIST3); $(RANLIB) meschach.a - -$(ZLIST1): $(HBASDE) zmatrix.h -zpart1: $(ZLIST1) - ar ru meschach.a $(ZLIST1); ranlib meschach.a - -$(ZLIST2): $(HBASE) zmatrix.h zmatrix2.h -zpart2: $(ZLIST2) - ar ru meschach.a $(ZLIST2); ranlib meschach.a - -$(OLDLIST): $(HBASE) sparse.h sparse2.h -oldpart: $(OLDLIST) - ar ru meschach.a $(OLDLIST); ranlib meschach.a - - - -####################################### - -tar: - - /bin/rm -f $(MES_PAK).tar - chmod 644 `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(OTHERS) $(HLIST) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` - chmod 755 configure - $(TAR) cvf $(MES_PAK).tar \ - `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(HLIST) $(OTHERS) \ - `echo $(TORTURE) | sed -e 's/\.o/.c/g'` \ - MACHINES DOC - -# use this only for PC machines -msdos-zip: - - /bin/rm -f $(MES_PAK).zip - chmod 644 `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(OTHERS) $(HLIST) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` - chmod 755 configure - $(ZIP) $(MES_PAK).zip \ - `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(HLIST) $(OTHERS) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` \ - MACHINES DOC - - -fullshar: - - /bin/rm -f $(MES_PAK).shar; - chmod 644 `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(OTHERS) $(HLIST) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` - chmod 755 configure - $(SHAR) `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(HLIST) $(OTHERS) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` \ - MACHINES DOC > $(MES_PAK).shar - -shar: - - /bin/rm -f meschach1.shar meschach2.shar meschach3.shar \ - meschach4.shar oldmeschach.shar meschach0.shar - chmod 644 `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(OTHERS) $(HLIST) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` - chmod 755 configure - $(SHAR) `echo $(LIST1) | sed -e 's/\.o/.c/g'` > meschach1.shar - $(SHAR) `echo $(LIST2) | sed -e 's/\.o/.c/g'` > meschach2.shar - $(SHAR) `echo $(LIST3) | sed -e 's/\.o/.c/g'` > meschach3.shar - $(SHAR) `echo $(ZLIST1) | sed -e 's/\.o/.c/g'` \ - `echo $(ZLIST2) | sed -e 's/\.o/.c/g'` > meschach4.shar - $(SHAR) `echo $(OLDLIST) | sed -e 's/\.o/.c/g'` > oldmeschach.shar - $(SHAR) $(OTHERS) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` \ - $(HLIST) DOC MACHINES > meschach0.shar - - -clean: - /bin/rm -f *.o core asx5213a.mat iotort.dat - -cleanup: - /bin/rm -f *.o core asx5213a.mat iotort.dat *.a - -alltorture: torture sptort ztorture memtort itertort mfuntort iotort - -torture:torture.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o torture torture.o \ - meschach.a $(LIBS) -sptort:sptort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o sptort sptort.o \ - meschach.a $(LIBS) -memtort: memtort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o memtort memtort.o \ - meschach.a $(LIBS) -ztorture:ztorture.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o ztorture ztorture.o \ - meschach.a $(LIBS) -itertort: itertort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o itertort itertort.o \ - meschach.a $(LIBS) - -iotort: iotort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o iotort iotort.o \ - meschach.a $(LIBS) -mfuntort: mfuntort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o mfuntort mfuntort.o \ - meschach.a $(LIBS) -tstmove: tstmove.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o tstmove tstmove.o \ - meschach.a $(LIBS) -tstpxvec: tstpxvec.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o tstpxvec tstpxvec.o \ - meschach.a $(LIBS) - //GO.SYSIN DD MACHINES/GCC/makefile echo MACHINES/GCC/machine.h 1>&2 sed >MACHINES/GCC/machine.h <<'//GO.SYSIN DD MACHINES/GCC/machine.h' 's/^-//' -/* machine.h. Generated automatically by configure. */ -/* Any machine specific stuff goes here */ -/* Add details necessary for your own installation here! */ - -/* This is for use with "configure" -- if you are not using configure - then use machine.van for the "vanilla" version of machine.h */ - -/* Note special macros: ANSI_C (ANSI C syntax) - SEGMENTED (segmented memory machine e.g. MS-DOS) - MALLOCDECL (declared if malloc() etc have - been declared) */ - - -#define ANSI_C 1 -#define NOT_SEGMENTED 1 -/* #undef HAVE_COMPLEX_H */ -#define HAVE_MALLOC_H 1 -#define STDC_HEADERS -#define HAVE_BCOPY 1 -#define HAVE_BZERO 1 -#define CHAR0ISDBL0 1 -#define WORDS_BIGENDIAN 1 -/* #undef U_INT_DEF */ - - -/* for basic or larger versions */ -#define COMPLEX 1 -#define SPARSE 1 - -/* for loop unrolling */ -/* #undef VUNROLL */ -/* #undef MUNROLL */ - -/* for segmented memory */ -#ifndef NOT_SEGMENTED -#define SEGMENTED -#endif - -/* if the system has malloc.h */ -#ifdef HAVE_MALLOC_H -#define MALLOCDECL 1 -#include -#endif - -/* any compiler should have this header */ -/* if not, change it */ -#include - - -/* Check for ANSI C memmove and memset */ -#ifdef STDC_HEADERS - -/* standard copy & zero functions */ -#define MEM_COPY(from,to,size) memmove((to),(from),(size)) -#define MEM_ZERO(where,size) memset((where),'\0',(size)) - -#ifndef ANSI_C -#define ANSI_C 1 -#endif - -#endif - -/* standard headers */ -#ifdef ANSI_C -#include -#include -#include -#include -#endif - - -/* if have bcopy & bzero and no alternatives yet known, use them */ -#ifdef HAVE_BCOPY -#ifndef MEM_COPY -/* nonstandard copy function */ -#define MEM_COPY(from,to,size) bcopy((char *)(from),(char *)(to),(int)(size)) -#endif -#endif - -#ifdef HAVE_BZERO -#ifndef MEM_ZERO -/* nonstandard zero function */ -#define MEM_ZERO(where,size) bzero((char *)(where),(int)(size)) -#endif -#endif - -/* if the system has complex.h */ -#ifdef HAVE_COMPLEX_H -#include -#endif - -/* If prototypes are available & ANSI_C not yet defined, then define it, - but don't include any header files as the proper ANSI C headers - aren't here */ -/* #undef HAVE_PROTOTYPES */ -#ifdef HAVE_PROTOTYPES -#ifndef ANSI_C -#define ANSI_C 1 -#endif -#endif - -/* 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 */ -#define REAL_DBL 1 - -/* 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 */ - -/******************** -#ifdef DBL_EPSILON -#define MACHEPS DBL_EPSILON -#endif -#ifdef M_MACHEPS -#ifndef MACHEPS -#define MACHEPS M_MACHEPS -#endif -#endif -********************/ - -#define M_MAX_INT 2147483647 -#ifdef M_MAX_INT -#ifndef MAX_RAND -#define MAX_RAND ((double)(M_MAX_INT)) -#endif -#endif - -/* for non-ANSI systems */ -#ifndef HUGE_VAL -#define HUGE_VAL HUGE -#endif - - -#ifdef ANSI_C -extern int isatty(int); -#endif - //GO.SYSIN DD MACHINES/GCC/machine.h mkdir MACHINES/RS6000 echo MACHINES/RS6000/machine.c 1>&2 sed >MACHINES/RS6000/machine.c <<'//GO.SYSIN DD MACHINES/RS6000/machine.c' 's/^-//' - -/************************************************************************** -** -** 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 matrix.a etc. - These are the routines that should be modified in order to take - full advantage of specialised architectures (pipelining, vector - processors etc). - */ -static char *rcsid = "$Header: /net/phoebus/epicsmgr/cvsroot/epics/extensions/src/SDDS/meschach/sourceFiles/m4,v 1.1.1.1 1999/04/14 14:16:22 borland Exp $"; - -#include "machine.h" - -/* __ip__ -- inner product */ -double __ip__(dp1,dp2,len) -register double *dp1, *dp2; -int len; -{ - register int len4; - register int i; - register double sum0, sum1, sum2, sum3; - - sum0 = sum1 = sum2 = sum3 = 0.0; - - len4 = len / 4; - len = len % 4; - - for ( i = 0; i < len4; i++ ) - { - sum0 += 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]; - } - sum0 += sum1 + sum2 + sum3; - dp1 += 4*len4; dp2 += 4*len4; - - for ( i = 0; i < len; i++ ) - sum0 += (*dp1++)*(*dp2++); - - return sum0; -} - -/* __mltadd__ -- scalar multiply and add c.f. v_mltadd() */ -void __mltadd__(dp1,dp2,s,len) -register double *dp1, *dp2, s; -register int len; -{ - register int i, 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; - - for ( i = 0; i < len; i++ ) - (*dp1++) += s*(*dp2++); -} - -/* __smlt__ scalar multiply array c.f. sv_mlt() */ -void __smlt__(dp,s,out,len) -register double *dp, s, *out; -register int len; -{ - register int i; - for ( i = 0; i < len; i++ ) - (*out++) = s*(*dp++); -} - -/* __add__ -- add arrays c.f. v_add() */ -void __add__(dp1,dp2,out,len) -register double *dp1, *dp2, *out; -register int len; -{ - register int i; - for ( i = 0; i < len; i++ ) - (*out++) = (*dp1++) + (*dp2++); -} - -/* __sub__ -- subtract arrays c.f. v_sub() */ -void __sub__(dp1,dp2,out,len) -register double *dp1, *dp2, *out; -register int len; -{ - register int i; - for ( i = 0; i < len; i++ ) - (*out++) = (*dp1++) - (*dp2++); -} - -/* __zero__ -- zeros an array of double precision numbers */ -void __zero__(dp,len) -register double *dp; -register int len; -{ - /* if a double precision zero is equivalent to a string of nulls */ - MEM_ZERO((char *)dp,len*sizeof(double)); - /* else, need to zero the array entry by entry */ - /************************************************* - while ( len-- ) - *dp++ = 0.0; - *************************************************/ -} - -/*********************************************************************** - ****** Faster versions ******** - ***********************************************************************/ - -/* __ip4__ -- compute 4 inner products in one go */ -void __ip4__(v0,v1,v2,v3,w,out,len) -double *v0, *v1, *v2, *v3, *w; -double out[4]; -int len; -{ - register int i, len2; - register double sum00, sum10, sum20, sum30, w_val0; - register double sum01, sum11, sum21, sum31, w_val1; - - len2 = len / 2; - len = len % 2; - sum00 = sum10 = sum20 = sum30 = 0.0; - sum01 = sum11 = sum21 = sum31 = 0.0; - for ( i = 0; i < len2; i++ ) - { - w_val0 = w[2*i]; - w_val1 = w[2*i+1]; - sum00 += v0[2*i] *w_val0; - sum01 += v0[2*i+1]*w_val1; - sum10 += v1[2*i] *w_val0; - sum11 += v1[2*i+1]*w_val1; - sum20 += v2[2*i] *w_val0; - sum21 += v2[2*i+1]*w_val1; - sum30 += v3[2*i] *w_val0; - sum31 += v3[2*i+1]*w_val1; - } - w += 2*len2; - v0 += 2*len2; - v1 += 2*len2; - v2 += 2*len2; - v3 += 2*len2; - for ( i = 0; i < len; i++ ) - { - w_val0 = w[i]; - sum00 += v0[i]*w_val0; - sum10 += v1[i]*w_val0; - sum20 += v2[i]*w_val0; - sum30 += v3[i]*w_val0; - } - out[0] = sum00 + sum01; - out[1] = sum10 + sum11; - out[2] = sum20 + sum21; - out[3] = sum30 + sum31; -} - -/* __lc4__ -- linear combinations: w <- w+a[0]*v0+ ... + a[3]*v3 */ -void __lc4__(v0,v1,v2,v3,w,a,len) -double *v0, *v1, *v2, *v3, *w; -double a[4]; -int len; -{ - register int i, len2; - register double a0, a1, a2, a3, tmp0, tmp1; - - len2 = len / 2; - len = len % 2; - - a0 = a[0]; a1 = a[1]; - a2 = a[2]; a3 = a[3]; - for ( i = 0; i < len2; i++ ) - { - tmp0 = w[2*i] + a0*v0[2*i]; - tmp1 = w[2*i+1] + a0*v0[2*i+1]; - tmp0 += a1*v1[2*i]; - tmp1 += a1*v1[2*i+1]; - tmp0 += a2*v2[2*i]; - tmp1 += a2*v2[2*i+1]; - tmp0 += a3*v3[2*i]; - tmp1 += a3*v3[2*i+1]; - w[2*i] = tmp0; - w[2*i+1] = tmp1; - } - w += 2*len2; - v0 += 2*len2; - v1 += 2*len2; - v2 += 2*len2; - v3 += 2*len2; - for ( i = 0; i < len; i++ ) - w[i] += a0*v0[i] + a1*v1[i] + a2*v2[i] + a3*v3[i]; -} - -/* __ma4__ -- multiply and add with 4 vectors: vi <- vi + ai*w */ -void __ma4__(v0,v1,v2,v3,w,a,len) -double *v0, *v1, *v2, *v3, *w; -double a[4]; -int len; -{ - register int i; - register double a0, a1, a2, a3, w0, w1, w2, w3; - - a0 = a[0]; a1 = a[1]; - a2 = a[2]; a3 = a[3]; - for ( i = 0; i < len; i++ ) - { - w0 = w[i]; - v0[i] += a0*w0; - v1[i] += a1*w0; - v2[i] += a2*w0; - v3[i] += a3*w0; - } -} //GO.SYSIN DD MACHINES/RS6000/machine.c echo MACHINES/RS6000/machine.h 1>&2 sed >MACHINES/RS6000/machine.h <<'//GO.SYSIN DD MACHINES/RS6000/machine.h' 's/^-//' - -/* Note special macros: ANSI_C (ANSI C syntax) - SEGMENTED (segmented memory machine e.g. MS-DOS) - MALLOCDECL (declared if malloc() etc have - been declared) */ - -#define ANSI_C 1 - -/* #undef MALLOCDECL */ -#define NOT_SEGMENTED 1 -/* #undef HAVE_COMPLEX_H */ -#define HAVE_MALLOC_H 1 -#define STDC_HEADERS 1 -#define HAVE_BCOPY 1 -#define HAVE_BZERO 1 -#define CHAR0ISDBL0 1 -#define WORDS_BIGENDIAN 1 -#define U_INT_DEF 1 - - -/* for basic or larger versions */ -#define COMPLEX 1 -#define SPARSE 1 - -/* for loop unrolling */ -/* #undef VUNROLL */ -/* #undef MUNROLL */ - -/* for segmented memory */ -#ifndef NOT_SEGMENTED -#define SEGMENTED -#endif - -/* if the system has malloc.h */ -#ifdef HAVE_MALLOC_H -#define MALLOCDECL 1 -#include -#endif - -/* any compiler should have this header */ -/* if not, change it */ -#include - - -/* Check for ANSI C memmove and memset */ -#ifdef STDC_HEADERS - -/* standard copy & zero functions */ -#define MEM_COPY(from,to,size) memmove((to),(from),(size)) -#define MEM_ZERO(where,size) memset((where),'\0',(size)) - -#ifndef ANSI_C -#define ANSI_C 1 -#endif - -#endif - -/* standard headers */ -#ifdef ANSI_C -#include -#include -#include -#include -#endif - - -/* if have bcopy & bzero and no alternatives yet known, use them */ -#ifdef HAVE_BCOPY -#ifndef MEM_COPY -/* nonstandard copy function */ -#define MEM_COPY(from,to,size) bcopy((char *)(from),(char *)(to),(int)(size)) -#endif -#endif - -#ifdef HAVE_BZERO -#ifndef MEM_ZERO -/* nonstandard zero function */ -#define MEM_ZERO(where,size) bzero((char *)(where),(int)(size)) -#endif -#endif - -/* if the system has complex.h */ -#ifdef HAVE_COMPLEX_H -#include -#endif - -/* If prototypes are available & ANSI_C not yet defined, then define it, - but don't include any header files as the proper ANSI C headers - aren't here */ -#define HAVE_PROTOTYPES 1 -#ifdef HAVE_PROTOTYPES -#ifndef ANSI_C -#define ANSI_C 1 -#endif -#endif - -/* 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 */ - -/******************** -#ifdef DBL_EPSILON -#define MACHEPS DBL_EPSILON -#endif -#ifdef M_MACHEPS -#ifndef MACHEPS -#define MACHEPS M_MACHEPS -#endif -#endif -********************/ - -#define M_MAX_INT 2147483647 -#ifdef M_MAX_INT -#ifndef MAX_RAND -#define MAX_RAND ((double)(M_MAX_INT)) -#endif -#endif - -/* for non-ANSI systems */ -#ifndef HUGE_VAL -#define HUGE_VAL HUGE -#endif - - -#ifdef ANSI_C -extern int isatty(int); -#endif - //GO.SYSIN DD MACHINES/RS6000/machine.h echo MACHINES/RS6000/makefile 1>&2 sed >MACHINES/RS6000/makefile <<'//GO.SYSIN DD MACHINES/RS6000/makefile' 's/^-//' -# Generated automatically from makefile.in by configure. -# -# Makefile for Meschach via autoconf -# -# Copyright (C) David Stewart & Zbigniew Leyk 1993 -# -# $Id: m4,v 1.1.1.1 1999/04/14 14:16:22 borland Exp $ -# - -srcdir = . -VPATH = . - -CC = cc - -DEFS = -DHAVE_CONFIG_H -LIBS = -lm -RANLIB = ranlib - - -CFLAGS = -O - - -.c.o: - $(CC) -c $(CFLAGS) $(DEFS) $< - -SHELL = /bin/sh -MES_PAK = mesch12a -TAR = tar -SHAR = stree -u -ZIP = zip -r -l -FLIST = FILELIST - -############################### - -LIST1 = copy.o err.o matrixio.o memory.o vecop.o matop.o pxop.o \ - submat.o init.o otherio.o machine.o matlab.o ivecop.o version.o \ - meminfo.o memstat.o -LIST2 = lufactor.o bkpfacto.o chfactor.o qrfactor.o solve.o hsehldr.o \ - givens.o update.o norm.o hessen.o symmeig.o schur.o svd.o fft.o \ - mfunc.o bdfactor.o -LIST3 = sparse.o sprow.o sparseio.o spchfctr.o splufctr.o \ - spbkp.o spswap.o iter0.o itersym.o iternsym.o -ZLIST1 = zmachine.o zcopy.o zmatio.o zmemory.o zvecop.o zmatop.o znorm.o \ - zfunc.o -ZLIST2 = zlufctr.o zsolve.o zmatlab.o zhsehldr.o zqrfctr.o \ - zgivens.o zhessen.o zschur.o - -# they are no longer supported -# if you use them add oldpart to all and sparse -OLDLIST = conjgrad.o lanczos.o arnoldi.o - -ALL_LISTS = $(LIST1) $(LIST2) $(LIST3) $(ZLIST1) $(ZLIST2) $(OLDLIST) - -HBASE = err.h meminfo.h machine.h matrix.h - -HLIST = $(HBASE) iter.h matlab.h matrix2.h oldnames.h sparse.h \ - sparse2.h zmatrix.h zmatrix2.h - -TORTURE = torture.o sptort.o ztorture.o memtort.o itertort.o \ - mfuntort.o iotort.o - -OTHERS = dmacheps.c extras.c fmacheps.c maxint.c makefile.in \ - README configure configure.in machine.h.in copyright \ - tutorial.c tutadv.c rk4.dat ls.dat makefile $(FLIST) - - -# Different configurations -all: part1 part2 part3 zpart1 zpart2 -basic: part1 part2 -sparse: part1 part2 part3 -complex: part1 part2 zpart1 zpart2 - - -$(LIST1): $(HBASE) -part1: $(LIST1) - ar ru meschach.a $(LIST1); $(RANLIB) meschach.a - -$(LIST2): $(HBASE) matrix2.h -part2: $(LIST2) - ar ru meschach.a $(LIST2); $(RANLIB) meschach.a -schur.o: schur.c $(HBASE) matrix2.h - cc -c $(DEFS) schur.c - -$(LIST3): $(HBASE) sparse.h sparse2.h -part3: $(LIST3) - ar ru meschach.a $(LIST3); $(RANLIB) meschach.a - -$(ZLIST1): $(HBASDE) zmatrix.h -zpart1: $(ZLIST1) - ar ru meschach.a $(ZLIST1); ranlib meschach.a - -$(ZLIST2): $(HBASE) zmatrix.h zmatrix2.h -zpart2: $(ZLIST2) - ar ru meschach.a $(ZLIST2); ranlib meschach.a - -$(OLDLIST): $(HBASE) sparse.h sparse2.h -oldpart: $(OLDLIST) - ar ru meschach.a $(OLDLIST); ranlib meschach.a - - - -####################################### - -tar: - - /bin/rm -f $(MES_PAK).tar - chmod 644 `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(OTHERS) $(HLIST) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` - chmod 755 configure - $(MAKE) list - $(TAR) cvf $(MES_PAK).tar \ - `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(HLIST) $(OTHERS) \ - `echo $(TORTURE) | sed -e 's/\.o/.c/g'` \ - MACHINES DOC - -# use this only for PC machines -msdos-zip: - - /bin/rm -f $(MES_PAK).zip - chmod 644 `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(OTHERS) $(HLIST) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` - chmod 755 configure - $(MAKE) list - $(ZIP) $(MES_PAK).zip \ - `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(HLIST) $(OTHERS) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` \ - MACHINES DOC - - -fullshar: - - /bin/rm -f $(MES_PAK).shar; - chmod 644 `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(OTHERS) $(HLIST) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` - chmod 755 configure - $(MAKE) list - $(SHAR) `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(HLIST) $(OTHERS) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` \ - MACHINES DOC > $(MES_PAK).shar - -shar: - - /bin/rm -f meschach1.shar meschach2.shar meschach3.shar \ - meschach4.shar oldmeschach.shar meschach0.shar - chmod 644 `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(OTHERS) $(HLIST) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` - chmod 755 configure - $(MAKE) list - $(SHAR) `echo $(LIST1) | sed -e 's/\.o/.c/g'` > meschach1.shar - $(SHAR) `echo $(LIST2) | sed -e 's/\.o/.c/g'` > meschach2.shar - $(SHAR) `echo $(LIST3) | sed -e 's/\.o/.c/g'` > meschach3.shar - $(SHAR) `echo $(ZLIST1) | sed -e 's/\.o/.c/g'` \ - `echo $(ZLIST2) | sed -e 's/\.o/.c/g'` > meschach4.shar - $(SHAR) `echo $(OLDLIST) | sed -e 's/\.o/.c/g'` > oldmeschach.shar - $(SHAR) $(OTHERS) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` \ - $(HLIST) DOC MACHINES > meschach0.shar - -list: - /bin/rm -f $(FLIST) - ls -lR `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - `echo $(TORTURE) | sed -e 's/\.o/.c/g'` \ - $(HLIST) $(OTHERS) MACHINES DOC \ - |awk '/^$$/ {print};/^[-d]/ {printf("%s %s %10d %s %s %s %s\n", \ - $$1,$$2,$$5,$$6,$$7,$$8,$$9)}; /^[^-d]/ {print}' \ - > $(FLIST) - - - -clean: - /bin/rm -f *.o core asx5213a.mat iotort.dat - -cleanup: - /bin/rm -f *.o core asx5213a.mat iotort.dat *.a - -alltorture: torture sptort ztorture memtort itertort mfuntort iotort - -torture:torture.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o torture torture.o \ - meschach.a $(LIBS) -sptort:sptort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o sptort sptort.o \ - meschach.a $(LIBS) -memtort: memtort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o memtort memtort.o \ - meschach.a $(LIBS) -ztorture:ztorture.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o ztorture ztorture.o \ - meschach.a $(LIBS) -itertort: itertort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o itertort itertort.o \ - meschach.a $(LIBS) - -iotort: iotort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o iotort iotort.o \ - meschach.a $(LIBS) -mfuntort: mfuntort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o mfuntort mfuntort.o \ - meschach.a $(LIBS) -tstmove: tstmove.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o tstmove tstmove.o \ - meschach.a $(LIBS) -tstpxvec: tstpxvec.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o tstpxvec tstpxvec.o \ - meschach.a $(LIBS) - //GO.SYSIN DD MACHINES/RS6000/makefile mkdir MACHINES/SPARC echo MACHINES/SPARC/machine.h 1>&2 sed >MACHINES/SPARC/machine.h <<'//GO.SYSIN DD MACHINES/SPARC/machine.h' 's/^-//' - -/* Note special macros: ANSI_C (ANSI C syntax) - SEGMENTED (segmented memory machine e.g. MS-DOS) - MALLOCDECL (declared if malloc() etc have - been declared) */ - -#define const - -/* #undef MALLOCDECL */ -#define NOT_SEGMENTED 1 -/* #undef HAVE_COMPLEX_H */ -#define HAVE_MALLOC_H 1 -/* #undef STDC_HEADERS */ -#define HAVE_BCOPY 1 -#define HAVE_BZERO 1 -#define CHAR0ISDBL0 1 -#define WORDS_BIGENDIAN 1 -/* #undef U_INT_DEF */ -#define VARARGS 1 - - -/* for basic or larger versions */ -#define COMPLEX 1 -#define SPARSE 1 - -/* for loop unrolling */ -/* #undef VUNROLL */ -/* #undef MUNROLL */ - -/* for segmented memory */ -#ifndef NOT_SEGMENTED -#define SEGMENTED -#endif - -/* if the system has malloc.h */ -#ifdef HAVE_MALLOC_H -#define MALLOCDECL 1 -#include -#endif - -/* any compiler should have this header */ -/* if not, change it */ -#include - - -/* Check for ANSI C memmove and memset */ -#ifdef STDC_HEADERS - -/* standard copy & zero functions */ -#define MEM_COPY(from,to,size) memmove((to),(from),(size)) -#define MEM_ZERO(where,size) memset((where),'\0',(size)) - -#ifndef ANSI_C -#define ANSI_C 1 -#endif - -#endif - -/* standard headers */ -#ifdef ANSI_C -#include -#include -#include -#include -#endif - - -/* if have bcopy & bzero and no alternatives yet known, use them */ -#ifdef HAVE_BCOPY -#ifndef MEM_COPY -/* nonstandard copy function */ -#define MEM_COPY(from,to,size) bcopy((char *)(from),(char *)(to),(int)(size)) -#endif -#endif - -#ifdef HAVE_BZERO -#ifndef MEM_ZERO -/* nonstandard zero function */ -#define MEM_ZERO(where,size) bzero((char *)(where),(int)(size)) -#endif -#endif - -/* if the system has complex.h */ -#ifdef HAVE_COMPLEX_H -#include -#endif - -/* If prototypes are available & ANSI_C not yet defined, then define it, - but don't include any header files as the proper ANSI C headers - aren't here */ -/* #undef HAVE_PROTOTYPES */ -#ifdef HAVE_PROTOTYPES -#ifndef ANSI_C -#define ANSI_C 1 -#endif -#endif - -/* 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 */ -#define REAL_DBL 1 - -/* 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 */ - -/******************** -#ifdef DBL_EPSILON -#define MACHEPS DBL_EPSILON -#endif -#ifdef M_MACHEPS -#ifndef MACHEPS -#define MACHEPS M_MACHEPS -#endif -#endif -********************/ - -#define M_MAX_INT 2147483647 -#ifdef M_MAX_INT -#ifndef MAX_RAND -#define MAX_RAND ((double)(M_MAX_INT)) -#endif -#endif - -/* for non-ANSI systems */ -#ifndef HUGE_VAL -#define HUGE_VAL HUGE -#endif - - -#ifdef ANSI_C -extern int isatty(int); -#endif - //GO.SYSIN DD MACHINES/SPARC/machine.h echo MACHINES/SPARC/makefile 1>&2 sed >MACHINES/SPARC/makefile <<'//GO.SYSIN DD MACHINES/SPARC/makefile' 's/^-//' -# # -# Makefile for Meschach for SUN SPARC cc -# -# Copyright (C) David Stewart & Zbigniew Leyk 1993 -# -# $Id: m4,v 1.1.1.1 1999/04/14 14:16:22 borland Exp $ -# - -srcdir = . -VPATH = . - -CC = cc - -DEFS = -DHAVE_CONFIG_H -LIBS = -lm -RANLIB = ranlib - - -CFLAGS = -O - - -.c.o: - $(CC) -c $(CFLAGS) $(DEFS) $< - -SHELL = /bin/sh -MES_PAK = mesch12a -TAR = tar -SHAR = stree -u -ZIP = zip -r -l - -############################### - -LIST1 = copy.o err.o matrixio.o memory.o vecop.o matop.o pxop.o \ - submat.o init.o otherio.o machine.o matlab.o ivecop.o version.o \ - meminfo.o memstat.o -LIST2 = lufactor.o bkpfacto.o chfactor.o qrfactor.o solve.o hsehldr.o \ - givens.o update.o norm.o hessen.o symmeig.o schur.o svd.o fft.o \ - mfunc.o bdfactor.o -LIST3 = sparse.o sprow.o sparseio.o spchfctr.o splufctr.o \ - spbkp.o spswap.o iter0.o itersym.o iternsym.o -ZLIST1 = zmachine.o zcopy.o zmatio.o zmemory.o zvecop.o zmatop.o znorm.o \ - zfunc.o -ZLIST2 = zlufctr.o zsolve.o zmatlab.o zhsehldr.o zqrfctr.o \ - zgivens.o zhessen.o zschur.o - -# they are no longer supported -# if you use them add oldpart to all and sparse -OLDLIST = conjgrad.o lanczos.o arnoldi.o - -ALL_LISTS = $(LIST1) $(LIST2) $(LIST3) $(ZLIST1) $(ZLIST2) $(OLDLIST) - - -HLIST = err.h iter.h machine.h matlab.h matrix.h matrix2.h \ - meminfo.h oldnames.h sparse.h sparse2.h \ - zmatrix.h zmatrix2.h - -TORTURE = torture.o sptort.o ztorture.o memtort.o itertort.o \ - mfuntort.o iotort.o - -OTHERS = dmacheps.c extras.c fmacheps.c maxint.c makefile.in \ - README configure configure.in machine.h.in copyright \ - tutorial.c tutadv.c rk4.dat ls.dat makefile - - -# Different configurations -all: part1 part2 part3 zpart1 zpart2 -basic: part1 part2 -sparse: part1 part2 part3 -complex: part1 part2 zpart1 zpart2 - - -HBASE = err.h meminfo.h machine.h matrix.h - -$(LIST1): $(HBASE) -part1: $(LIST1) - ar ru meschach.a $(LIST1); $(RANLIB) meschach.a - -$(LIST2): $(HBASE) matrix2.h -part2: $(LIST2) - ar ru meschach.a $(LIST2); $(RANLIB) - -$(LIST3): $(HBASE) sparse.h sparse2.h -part3: $(LIST3) - ar ru meschach.a $(LIST3); $(RANLIB) meschach.a - -$(ZLIST1): $(HBASDE) zmatrix.h -zpart1: $(ZLIST1) - ar ru meschach.a $(ZLIST1); ranlib meschach.a - -$(ZLIST2): $(HBASE) zmatrix.h zmatrix2.h -zpart2: $(ZLIST2) - ar ru meschach.a $(ZLIST2); ranlib meschach.a - -$(OLDLIST): $(HBASE) sparse.h sparse2.h -oldpart: $(OLDLIST) - ar ru meschach.a $(OLDLIST); ranlib meschach.a - - - -####################################### - -tar: - - /bin/rm -f $(MES_PAK).tar - chmod 644 `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(OTHERS) $(HLIST) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` - chmod 755 configure - $(TAR) cvf $(MES_PAK).tar \ - `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(HLIST) $(OTHERS) \ - `echo $(TORTURE) | sed -e 's/\.o/.c/g'` \ - MACHINES DOC - -# use this only for PC machines -msdos-zip: - - /bin/rm -f $(MES_PAK).zip - chmod 644 `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(OTHERS) $(HLIST) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` - chmod 755 configure - $(ZIP) $(MES_PAK).zip \ - `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(HLIST) $(OTHERS) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` \ - MACHINES DOC - - -fullshar: - - /bin/rm -f $(MES_PAK).shar; - chmod 644 `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(OTHERS) $(HLIST) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` - chmod 755 configure - $(SHAR) `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(HLIST) $(OTHERS) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` \ - MACHINES DOC > $(MES_PAK).shar - -shar: - - /bin/rm -f meschach1.shar meschach2.shar meschach3.shar \ - meschach4.shar oldmeschach.shar meschach0.shar - chmod 644 `echo $(ALL_LISTS) | sed -e 's/\.o/.c/g'` \ - $(OTHERS) $(HLIST) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` - chmod 755 configure - $(SHAR) `echo $(LIST1) | sed -e 's/\.o/.c/g'` > meschach1.shar - $(SHAR) `echo $(LIST2) | sed -e 's/\.o/.c/g'` > meschach2.shar - $(SHAR) `echo $(LIST3) | sed -e 's/\.o/.c/g'` > meschach3.shar - $(SHAR) `echo $(ZLIST1) | sed -e 's/\.o/.c/g'` \ - `echo $(ZLIST2) | sed -e 's/\.o/.c/g'` > meschach4.shar - $(SHAR) `echo $(OLDLIST) | sed -e 's/\.o/.c/g'` > oldmeschach.shar - $(SHAR) $(OTHERS) `echo $(TORTURE) | sed -e 's/\.o/.c/g'` \ - $(HLIST) DOC MACHINES > meschach0.shar - - -clean: - /bin/rm -f *.o core asx5213a.mat iotort.dat - -cleanup: - /bin/rm -f *.o core asx5213a.mat iotort.dat *.a - -alltorture: torture sptort ztorture memtort itertort mfuntort iotort - -torture:torture.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o torture torture.o \ - meschach.a $(LIBS) -sptort:sptort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o sptort sptort.o \ - meschach.a $(LIBS) -memtort: memtort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o memtort memtort.o \ - meschach.a $(LIBS) -ztorture:ztorture.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o ztorture ztorture.o \ - meschach.a $(LIBS) -itertort: itertort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o itertort itertort.o \ - meschach.a $(LIBS) - -iotort: iotort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o iotort iotort.o \ - meschach.a $(LIBS) -mfuntort: mfuntort.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o mfuntort mfuntort.o \ - meschach.a $(LIBS) -tstmove: tstmove.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o tstmove tstmove.o \ - meschach.a $(LIBS) -tstpxvec: tstpxvec.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o tstpxvec tstpxvec.o \ - meschach.a $(LIBS) - //GO.SYSIN DD MACHINES/SPARC/makefile mkdir MACHINES/Linux echo MACHINES/Linux/makefile 1>&2 sed >MACHINES/Linux/makefile <<'//GO.SYSIN DD MACHINES/Linux/makefile' 's/^-//' -# Generated automatically from makefile.in by configure. -# -# Makefile for Meschach via autoconf -# -# Copyright (C) David Stewart & Zbigniew Leyk 1993 -# -# $Id: m4,v 1.1.1.1 1999/04/14 14:16:22 borland Exp $ -# - -srcdir = . -VPATH = . - -CC = cc - -DEFS = -DHAVE_CONFIG_H -LIBS = -lm -RANLIB = ranlib - - -CFLAGS = -O - - -.c.o: - $(CC) -c $(CFLAGS) $(DEFS) $< - -SHELL = /bin/sh -MES_PAK = mesch12a -TAR = tar -SHAR = stree -u -ZIP = zip -r -l -FLIST = FILELIST bigmail CUT HERE............ test -w meschach0.shar && test -r 24048P00 && test -r 24048P01 && test -r 24048P02 && test -r 24048P03 && ( cat 24048P00 >> meschach0.shar; rm 24048P00 cat 24048P01 >> meschach0.shar; rm 24048P01 cat 24048P02 >> meschach0.shar; rm 24048P02 cat 24048P03 >> meschach0.shar; rm 24048P03 )