From 6589eb7f7c11fe514b5a25b0dafe1d9672c0a8c9 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 18 Jan 2012 23:44:50 +0000 Subject: [PATCH] Removed some unused numerical integration utilities --- Cantera/src/base/mdp_allo.cpp | 1784 ---------------------- Cantera/src/base/mdp_allo.h | 734 --------- Cantera/src/numerics/NonlinearSolver.cpp | 963 ------------ Cantera/src/numerics/NonlinearSolver.h | 329 ---- Cantera/src/numerics/ResidEval.h | 102 -- Cantera/src/numerics/ResidJacEval.cpp | 282 ---- Cantera/src/numerics/ResidJacEval.h | 175 --- Cantera/src/numerics/SquareMatrix.cpp | 114 -- Cantera/src/numerics/SquareMatrix.h | 91 -- 9 files changed, 4574 deletions(-) delete mode 100644 Cantera/src/base/mdp_allo.cpp delete mode 100644 Cantera/src/base/mdp_allo.h delete mode 100644 Cantera/src/numerics/NonlinearSolver.cpp delete mode 100644 Cantera/src/numerics/NonlinearSolver.h delete mode 100755 Cantera/src/numerics/ResidEval.h delete mode 100644 Cantera/src/numerics/ResidJacEval.cpp delete mode 100644 Cantera/src/numerics/ResidJacEval.h delete mode 100644 Cantera/src/numerics/SquareMatrix.cpp delete mode 100644 Cantera/src/numerics/SquareMatrix.h diff --git a/Cantera/src/base/mdp_allo.cpp b/Cantera/src/base/mdp_allo.cpp deleted file mode 100644 index 5e0785dd1..000000000 --- a/Cantera/src/base/mdp_allo.cpp +++ /dev/null @@ -1,1784 +0,0 @@ -/** - * @file mdp_allo.cpp - * Definitions for dynamic allocation of multidimensional pointer arrays - */ -/* - * Copywrite 2004 Sandia Corporation. Under the terms of Contract - * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government - * retains certain rights in this software. - * See file License.txt for licensing information. - */ -#include - -#include -#include -#include -#include - -#include "mdp_allo.h" - -using namespace std; - -namespace mdp { - /* - * Allocate global storage for 2 debugging ints that are used in IO of - * error information. - */ -#ifdef MDP_MPDEBUGIO - int MDP_MP_Nprocs = 1; - int MDP_MP_myproc = 0; -#endif - /* - * Error Handling - * 7 print and exit - * 6 exit - * 5 print and create a divide by zero for stack trace analysis. - * 4 create a divide by zero for stack analysis trace - * 3 print a message and throw the bad_alloc exception. - * 2 throw the bad_alloc exception and be quite - * 1 print a message and return from package with the NULL pointer - * 0 Keep completely silent about the matter and return with - * a null pointer. - * - * -> Right now, the only way to change this option is to right here - */ - int MDP_ALLO_errorOption = 3; - - inline int MinI(const int &x, const int &y) { - return ( ( x < y ) ? x : y ); - } - inline int MaxI(const int &x, const int &y) { - return ( ( x > y ) ? x : y ); - } - - - const int MDP_ALLOC_INTERFACE_ERROR = -230346; - - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - static void mdp_alloc_eh(const char * const rname, const int bytes) - /************************************************************************* - * - * mdp_alloc_eh: - * - * Error Handling - * 7 print and exit - * 6 exit - * 5 print and create a divide by zero for stack trace analysis. - * 4 create a divide by zero for stack trace analysis. - * 3 print a message and throw the bad_alloc exception. - * 2 throw the bad_alloc exception and be quite - * 1 print a message and return from package with the NULL pointer - * 0 Keep completely silent about the matter and return with - * a null pointer. - **************************************************************************/ - { - double cd = 0.0; - static char mesg[64]; - if (bytes == MDP_ALLOC_INTERFACE_ERROR) { -#ifdef MDP_MPDEBUGIO - sprintf(mesg,"MDP_ALLOC Interface ERROR P_%d: %s", MDP_MP_my_proc, - rname); -#else - sprintf(mesg,"MDP_ALLOC Interface ERROR: %s", rname); -#endif - } else { - sprintf(mesg,"%s ERROR: out of memory while mallocing %d bytes", - rname, bytes); - } - if (MDP_ALLO_errorOption % 2 == 1) { - fprintf(stderr, "\n%s", mesg); -#ifdef MDP_MPDEBUGIO - if (MDP_MP_Nprocs > 1) { - fprintf(stderr,": proc = %d\n", MDP_MP_myproc); - } else { - fprintf(stderr,"\n"); - } -#else - fprintf(stderr,"\n"); -#endif - } - fflush(stderr); - if (MDP_ALLO_errorOption == 2 || MDP_ALLO_errorOption == 3) { - throw std::bad_alloc(); - } - if (MDP_ALLO_errorOption == 4 || MDP_ALLO_errorOption == 5) cd = 1.0 / cd; - if (MDP_ALLO_errorOption > 5) { - std::exit(-1); - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - static void mdp_alloc_eh2(const char * const rname) - - /************************************************************************* - * - * mdp_alloc_eh2: - * - * Second Level Error Handling - * This routine is used at the second level. - * It will be called within the routine whenever another routine - * returns a NULL pointer. - * - **************************************************************************/ - { - if (MDP_ALLO_errorOption == 1 || - MDP_ALLO_errorOption == 3 || - MDP_ALLO_errorOption == 5 || - MDP_ALLO_errorOption == 7) { -#ifdef MDP_MPDEBUGIO - if (MDP_MP_Nprocs > 1) { - fprintf(stderr,": proc = %d, %s ERROR: returning with null pointer\n", - MDP_MP_myproc, rname); - } else { - fprintf(stderr,"%s ERROR: returning with null pointer", rname); - } -#else - fprintf(stderr,"%s ERROR: returning with null pointer", rname); -#endif - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - /****************************************************************************/ -#ifndef HAVE_ARRAY_ALLOC - /****************************************************************************/ - static double *smalloc(size_t n); - - /***************************************************************************** - * - * Dynamic Allocation of Multidimensional Arrays - *---------------------------------------------------------------------------- - * - * Example Usage: - * - * typedef struct - * { int bus1; - * int bus2; - * int dest; - * } POINT; - * - * POINT **points, corner; - * - * points = (POINT **) array_alloc (2, x, y, sizeof(POINT)); - * ^ ^ ^ - * | | | - * number of dimensions----------+ | | - * | | - * first dimension max-------------+ | - * | - * second dimension max----------------+ - * - * (points may be now be used as if it were declared - * POINT points[x][y]) - * - * This particular version is limited to dimensions of 3 or less. - * - * corner = points[2][3]; (refer to the structure as you would any array) - * - * free (points); (frees the entire structure in one fell swoop) - * - ****************************************************************************/ - /***************************************************************************** - * The following section is a commented section containing - * an example main code: - ****************************************************************************** - *double *array_alloc(); - *main() - *{ - * int ***temp; - * int *temp2; - * int i, j, k; - * int il, jl, kl; - * - * malloc_debug(2); - * il = 2; - * jl = 3; - * kl = 3; - * temp = (int ***) array_alloc(3,il,jl,kl,sizeof(int)); - * for (i=0; i 4) { - (void) fprintf(stderr, - "mdp_array_alloc ERROR: number of dimensions, %d, is > 4\n", - numdim); - return NULL; - } - - dim[0].index = va_arg(va, int); - - if (dim[0].index <= 0) { -#ifdef MDP_MEMDEBUG - (void) fprintf(stderr, "WARNING: mdp_array_alloc called with first " - "dimension <= 0, %d\n\twill return the nil pointer\n", - (int) (dim[0].index)); -#endif - return((double *) NULL); - } - - dim[0].total = dim[0].index; - dim[0].size = sizeof(void *); - dim[0].off = 0; - for (i = 1; i < numdim; i++) { - dim[i].index = va_arg(va, int); - if (dim[i].index <= 0) { - (void) fprintf(stderr, - "WARNING: mdp_array_alloc called with dimension %d <= 0, " - "%d\n", i+1, (int) (dim[i].index)); - fprintf(stderr, "\twill return the nil pointer\n"); - return((double *) NULL); - } - dim[i].total = dim[i-1].total * dim[i].index; - dim[i].size = sizeof(void *); - dim[i].off = dim[i-1].off + dim[i-1].total * dim[i-1].size; - } - - dim[numdim-1].size = va_arg(va, int); - va_end(va); - - /* - * Round up the last offset value so data is properly aligned. - */ - - dim[numdim-1].off = dim[numdim-1].size * - ((dim[numdim-1].off+dim[numdim-1].size-1)/dim[numdim-1].size); - - total = dim[numdim-1].off + dim[numdim-1].total * dim[numdim-1].size; - - dfield = (double *) smalloc((size_t) total); - field = (char *) dfield; - - for (i = 0; i < numdim - 1; i++) { - ptr = (char **) (field + dim[i].off); - data = (char *) (field + dim[i+1].off); - for (j = 0; j < dim[i].total; j++) { - ptr[j] = data + j * dim[i+1].size * dim[i+1].index; - } - } - - return dfield; - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - static double *smalloc(size_t n) - - /************************************************************************** - * smalloc: safe version of malloc - * - * This version of smalloc assigns space in even chunks of 8 bytes only - * - **************************************************************************/ - { -#ifdef MDP_MEMDEBUG - static int firsttime = 1; - FILE *file; -#endif - double *pntr; - if (n < 0) { - if (MDP_ALLO_errorOption == 7 || - MDP_ALLO_errorOption == 5 || MDP_ALLO_errorOption == 3 || - MDP_ALLO_errorOption == 1) { - (void) fprintf(stderr, "smalloc ERROR: Non-positive argument. (%d)\n", (int) n); - return NULL; - } - } - else if (n == 0) pntr = NULL; - else { - n = ((n - 1) / 8); - n = (n + 1) * 8; - pntr = (double *) malloc((size_t) n); - } - if (pntr == NULL && n != 0) { - if (MDP_ALLO_errorOption == 7 || - MDP_ALLO_errorOption == 5 || MDP_ALLO_errorOption == 3 || - MDP_ALLO_errorOption == 1) { - fprintf(stderr, "smalloc : Out of space - number of bytes " - "requested = %d\n", (int) n); - } - } -#ifdef MDP_MEMDEBUG - if (firsttime) { - firsttime = 0; - file = fopen("memops.txt", "w"); - } else { - file = fopen("memops.txt", "a"); - } - fprintf(file, "%x %d malloc\n", pntr, n); - fclose(file); -#endif - return pntr; - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_safe_free (void **ptr) - - /************************************************************************* - * - * mdp_safe_free(): - * - * This version of free calls the system's free function - * with maximum error checking. It also doesn't call free if ptr is - * the NULL pointer already. - * It will then set the freed pointer to NULL. Thus, a convention may - * be established wherein all pointers that can be malloced can be - * set to NULL if they are not malloced. - **************************************************************************/ - { -#ifdef MDP_MEMDEBUG - FILE *file; -#endif - if (ptr == NULL) { - mdp_alloc_eh("mdp_safe_free: handle is NULL", MDP_ALLOC_INTERFACE_ERROR); - } - if (*ptr != NULL) { -#ifdef MDP_MEMDEBUG - file = fopen("memops.txt", "a"); - Fprintf(file, "%x free\n", *ptr); - fflush(file); - if ( (int) *ptr == 0x00000001) { - Fprintf(stderr, "FOUND IT!\n"); - exit(-1); - } - fclose(file); -#endif - free(*ptr); - /* - * Set the value of ptr to NULL, so that further references - * to it will be flagged. - */ - *ptr = NULL; - } - } - /****************************************************************************/ -#endif - /***************************************************************************** - * - * Wrapper Functions - * -------------------- - * - * The function definitions below are wrappers around array_alloc for - * common operations. The following principles are followed: - * - * Argument dimensions le 0 are increased to 1 before calling array_alloc. - * Thus, something is always malloced during a call. The reason for this is - * that it minimizes the number of special cases in the calling program. - * - * A pointer to something else other than NULL indicates that that pointer - * has been previously malloced. Thus, it can be freed. Note, after a free - * operation, this package always sets the pointer to NULL before returning. - * - * "safe_alloc" routines try to free the pointer if nonNULL, before calling - * the base alloc_int_#() routines. - * - * The regular routines initialize the malloced space, unless told to not - * do so. - * - * The function memcpy and memset are used where possible to increase speed. - * - * Naming Convention: - * - * The names of functions have the following format: - * - * mdp_[safe]_[operation]_[underlyingDataType]_[Dimension]() - * - * safe -> optional name added to the function name. These - * functions differ in that a handle to the pointer is passed - * to the function instead of just the pointer being - * returned. A check to see whether the pointer has been - * previously malloced is made, before a new allocation - * is made. If previously malloced, the pointer is freed - * first. - * - * operation -> name of the operation, such as alloc, realloc, - * init, or copy. - * - * underlyingDataType -> dble, int, ptr, or char - * - * dimension -> 1 or 2. For dimensions greater than 1, the data - * is laid out in F77 or blas-compatible format - * where data is contiguous wrt rows being the - * inner loop. - * - *****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - int *mdp_alloc_int_1(int nvalues, const int val) - - /************************************************************************** - * - * mdp_alloc_int_1: - * - * Allocate and initialize a one dimensional array of integers. - * - * Input - * ------- - * nvalues = Length of the array - * val = intialization value - * Return - * ------ - * Pointer to the intialized integer array - * Failures are indicated by returning the NULL pointer. - **************************************************************************/ - { - int *array; - if (nvalues <= 0) nvalues = 1; - array= (int *) mdp_array_alloc(1, nvalues, sizeof(int)); - if (array != NULL) { - if (val != MDP_INT_NOINIT) { - if (val == 0) { - (void) memset(array, 0, sizeof(int)*nvalues); - } else { - for (int i = 0; i < nvalues; i++) array[i] = val; - } - } - } else { - mdp_alloc_eh("mdp_alloc_int_1", nvalues * sizeof(int)); - } - return array; - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_safe_alloc_int_1(int **array_hdl, int nvalues, const int val) - - /************************************************************************* - * - * mdp_safe_alloc_int_1: - * - * Allocates and/or initializse a one dimensional array of integers. - * - * Input - * ------- - * *array_hdl = Previous value of pointer. If non-NULL will try - * to free the memory at this address. - * nvalues = Length of the array - * val = intialization value - * Output - * ------ - * *array_hdl = This value is initialized to the correct address - * of the array. - * A NULL value in the position indicates an error. - **************************************************************************/ - { - if (array_hdl == NULL) { - mdp_alloc_eh("mdp_safe_alloc_int_1: handle is NULL", - MDP_ALLOC_INTERFACE_ERROR); - return; - } - if (*array_hdl != NULL) mdp_safe_free((void **) array_hdl); - *array_hdl = mdp_alloc_int_1(nvalues, val); - if (*array_hdl == NULL) mdp_alloc_eh2("mdp_safe_alloc_int_1"); - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void - mdp_realloc_int_1(int **array_hdl, int new_length, int old_length, - const int defval) - - /************************************************************************* - * - * mdp_realloc_int_1_(array_hdl, new_num_ptrs, old_num_ptrs); - * - * Reallocates a one dimensional array of ints. - * This routine always allocates space for at least one int. - * Calls the smalloc() routine to ensure that all malloc - * calls go through one location. This routine will then copy - * the pertinent information from the old array to the - * new array. - * - * Input - * ------- - * array_hdl = Pointer to the global variable that - * holds the old and (eventually new) - * address of the array of integers to be reallocated - * new_length = Length of the array - * old_length = Length of the old array - **************************************************************************/ - { - if (new_length == old_length) return; - if (new_length <= 0) { -#ifdef MDP_MPDEBUGIO - fprintf(stderr, - "Warning: mdp_realloc_int_1 P_%d: called with n = %d\n", - MDP_MP_myproc, new_length); -#else - fprintf(stderr, - "Warning: mdp_realloc_int_1: called with n = %d\n", - new_length); -#endif - new_length = 1; - } - if (old_length < 0) old_length = 0; - if (new_length == old_length) return; - size_t bytenum = new_length * sizeof(int); - int *array = (int *) smalloc(bytenum); - if (array != NULL) { - if (*array_hdl) { - if (old_length > 0) bytenum = sizeof(int) * old_length; - else bytenum = 0; - if (new_length < old_length) bytenum = sizeof(int) * new_length; - if (bytenum > 0) { - (void) memcpy((void *) array, (const void *) *array_hdl, bytenum); - } - mdp_safe_free((void **) array_hdl); - } else { - old_length = 0; - } - *array_hdl = array; - if ((defval != MDP_INT_NOINIT) && (new_length > old_length)) { - if (defval == 0) { - bytenum = sizeof(int) * (new_length - old_length); - (void) memset((void *)(array+old_length), 0, bytenum); - } else { - for (int i = old_length; i < new_length; i++) { - array[i] = defval; - } - } - } - } else { - mdp_alloc_eh("mdp_realloc_int_1", static_cast(bytenum)); - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - int **mdp_alloc_int_2(int ndim1, int ndim2, const int val) - - /************************************************************************* - * - * mdp_alloc_int_2: - * - * Allocate and initialize a two dimensional array of ints. - * - * Input - * ------- - * ndim1 = Length of the first dimension of the array - * ndim2 = Length of the second dimension of the array - * val = intialization value - * Return - * ------ - * Pointer to the intialized integer array - * Failures are indicated by returning the NULL pointer. - **************************************************************************/ - { - int i; - int **array, *dptr; - if (ndim1 <= 0) ndim1 = 1; - if (ndim2 <= 0) ndim2 = 1; - array = (int **) mdp_array_alloc(2, ndim1, ndim2, sizeof(int)); - if (array != NULL) { - if (val != MDP_INT_NOINIT) { - if (val == 0) { - (void) memset((void *) array[0], 0, ndim1 * ndim2 * sizeof(int)); - } else { - dptr = &(array[0][0]); - for (i = 0; i < ndim1 * ndim2; i++) dptr[i] = val; - } - } - } else { - mdp_alloc_eh("mdp_alloc_int_2", - sizeof(int) * ndim1 * ndim2 + - ndim1 * sizeof(void *)); - } - return array; - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - double *mdp_alloc_dbl_1(int nvalues, const double val) - - /************************************************************************* - * - * mdp_alloc_dbl_1: - * - * Allocate and initialize a one dimensional array of doubles. - * - * Input - * ------- - * nvalues = Length of the array - * val = intialization value - * Return - * ------ - * Pointer to the intialized array of doubles - * Failures are indicated by returning the NULL pointer. - **************************************************************************/ - { - int i; - double *array; - if (nvalues <= 0) nvalues = 1; - array = (double *) mdp_array_alloc(1, nvalues, sizeof(double)); - if (array != NULL) { - if (val != MDP_DBL_NOINIT) { - if (val == 0.0) { - (void) memset((void *) array, 0, nvalues * sizeof(double)); - } else { - for (i = 0; i < nvalues; i++) array[i] = val; - } - } - } else { - mdp_alloc_eh("mdp_alloc_dbl_1", nvalues * sizeof(double)); - } - return array; - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_safe_alloc_dbl_1(double **array_hdl, int nvalues, const double val) - - /************************************************************************* - * - * mdp_safe_alloc_dbl_1: - * - * Allocates and/or initializse a one dimensional array of doubles. - * - * Input - * ------- - * *array_hdl = Previous value of pointer. If non-NULL will try - * to free the memory at this address. - * nvalues = Length of the array - * val = intialization value - * Output - * ------ - * *array_hdl = This value is initialized to the correct address - * of the array. - * A NULL value in the position indicates an error. - **************************************************************************/ - { - if (array_hdl == NULL) { - mdp_alloc_eh("mdp_safe_alloc_dbl_1: handle is NULL", - MDP_ALLOC_INTERFACE_ERROR); - return; - } - if (*array_hdl != NULL) mdp_safe_free((void **) array_hdl); - *array_hdl = mdp_alloc_dbl_1(nvalues, val); - if (*array_hdl == NULL) mdp_alloc_eh2("mdp_safe_alloc_dbl_1"); - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_realloc_dbl_1(double **array_hdl, int new_length, - int old_length, const double defval) - - /************************************************************************* - * - * mdp_realloc_dbl_1_(array_hdl, new_num_ptrs, old_num_ptrs); - * - * Reallocates a one dimensional array of doubles. - * This routine always allocates space for at least one dbl. - * Calls the smalloc() routine to ensure that all malloc - * calls go through one location. This routine will then copy - * the pertinent information from the old array to the - * new array. - * - * Input - * ------- - * array_hdl = Pointer to the global variable that - * holds the old and (eventually new) - * address of the array of doubles to be reallocated - * new_length = Length of the array - * old_length = Length of the old array - **************************************************************************/ - { - if (new_length == old_length) return; - if (new_length <= 0) { -#ifdef MDP_MPDEBUGIO - fprintf(stderr, "Warning: mdp_realloc_dbl_1 P_%d: called with n = %d\n", - MDP_MP_myproc, new_length); -#else - fprintf(stderr, "Warning: mdp_realloc_dbl_1: called with n = %d\n", - new_length); -#endif - new_length = 1; - } - if (old_length < 0) old_length = 0; - if (new_length == old_length) return; - size_t bytenum = new_length * sizeof(double); - double *array = (double *) smalloc(bytenum); - if (array != NULL) { - if (*array_hdl) { - if (old_length > 0) bytenum = sizeof(double) * old_length; - else bytenum = 0; - if (new_length < old_length) bytenum = sizeof(double) * new_length; - (void) memcpy((void *) array, (const void *) *array_hdl, bytenum); - mdp_safe_free((void **) array_hdl); - } else { - old_length = 0; - } - *array_hdl = array; - if ((defval != MDP_DBL_NOINIT) && (new_length > old_length)) { - if (defval == 0) { - bytenum = sizeof(double) * (new_length - old_length); - (void) memset((void *)(array+old_length), 0, bytenum); - } else { - for (int i = old_length; i < new_length; i++) { - array[i] = defval; - } - } - } - } else { - mdp_alloc_eh("mdp_realloc_dbl_1", static_cast(bytenum)); - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - double **mdp_alloc_dbl_2(int ndim1, int ndim2, const double val) - - /************************************************************************* - * - * mdp_alloc_dbl_2: - * - * Allocate a two dimensional array of doubles. The array is in - * fortran order and can be accessed via the following form: - * - * dblArray[ndim1][ndim2] - * - * Note, ndim2 is the inner dimension. i.e., the array is - * in column ordering. - * - * Input - * ------- - * ndim1 = Length of the first dimension of the array - * ndim2 = Length of the second dimension of the array - * val = intialization value - * Return - * ------ - * Pointer to the intialized integer array - * Failures are indicated by returning the NULL pointer. - **************************************************************************/ - { - int i; - double **array, *dptr; - if (ndim1 <= 0) ndim1 = 1; - if (ndim2 <= 0) ndim2 = 1; - array = (double **) mdp_array_alloc(2, ndim1, ndim2, sizeof(double)); - if (array != NULL) { - if (val != MDP_DBL_NOINIT) { - if (val == 0.0) { - (void) memset((void *) array[0], 0, ndim1*ndim2 * sizeof(double)); - } else { - dptr = &(array[0][0]); - for (i = 0; i < ndim1*ndim2; i++) dptr[i] = val; - } - } - } else { - mdp_alloc_eh("mdp_alloc_dbl_2", - sizeof(double) * ndim1 * ndim2 + - ndim1 * sizeof(void *)); - } - return array; - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - char *mdp_alloc_char_1(int nvalues, const char val) - - /************************************************************************* - * - * mdp_alloc_char_1: - * - * Allocate and initialize a one dimensional array of characters. - * - * Input - * ------- - * nvalues = Length of the array - * val = intialization value - * Return - * ------ - * Pointer to the intialized character array - * Failures are indicated by returning the NULL pointer. - **************************************************************************/ - { - int i; - char *array; - if (nvalues <= 0) nvalues = 1; - array = (char *) mdp_array_alloc(1, nvalues, sizeof(char)); - if (array != NULL) { - if (val == '\0') { - (void) memset((void *)array, 0, sizeof(char)*nvalues); - } else { - for (i = 0; i < nvalues; i++) array[i] = val; - } - } else { - mdp_alloc_eh("mdp_alloc_char_1", nvalues * sizeof(char)); - } - return array; - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_safe_alloc_char_1(char **array_hdl, int nvalues, const char val) - - /************************************************************************* - * - * mdp_safe_alloc_char_1: - * - * Allocate and initialize a one dimensional array of characters, - * deallocating the space before hand. - * - * This routine will free any old memory that was located at that - * position, before it will allocate a new vector. - * Allocates and/or initializse a one dimensional array of characters. - * - * Input - * ------- - * array_hdl = Previous value of pointer. If non-NULL will try - * to free the memory at this address. - * nvalues = Length of the array - * val = intialization value - * Output - * ------ - * *array_hdl = This value is initialized to the correct address - * of the array. - * A NULL value in the position indicates an error. - **************************************************************************/ - { - if (array_hdl == NULL) { - mdp_alloc_eh("mdp_safe_alloc_char_1: handle is NULL", - MDP_ALLOC_INTERFACE_ERROR); - return; - } - if (*array_hdl != NULL) mdp_safe_free((void **) array_hdl); - *array_hdl = mdp_alloc_char_1(nvalues, val); - if (*array_hdl == NULL) mdp_alloc_eh2("mdp_safe_alloc_char_1"); - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - - - void mdp_safe_alloc_dbl_2(double ***array_hdl, int ndim1, int ndim2, - const double val) - - /************************************************************************* - * - * mdp_safe_alloc_dbl_2: - * - * Allocate and initialize a two dimensional array of doubles. - * - * Input - * ------- - * *array_hdl = Previous value of pointer. If non-NULL will try - * to free the memory at this address. - * ndim1 = Length of the array - * ndim2 = Length of inner loop of the array - * val = intialization value - * Return - * ------ - * *array_hdl = This value is initialized to the correct address - * of the array. - * A NULL value in the position indicates an error. - **************************************************************************/ - { - if (array_hdl == NULL) { - mdp_alloc_eh("mdp_safe_alloc_dbl_2: handle is NULL", - MDP_ALLOC_INTERFACE_ERROR); - return; - } - if (*array_hdl != NULL) mdp_safe_free((void **) array_hdl); - *array_hdl = mdp_alloc_dbl_2(ndim1, ndim2, val); - if (*array_hdl == NULL) mdp_alloc_eh2("mdp_safe_alloc_dbl_2"); - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_realloc_dbl_2(double ***array_hdl, int ndim1, int ndim2, - int ndim1Old, int ndim2Old, const double val) - - /************************************************************************* - * - * mdp_realloc_dbl_2: - * - * mdp_realloc_dbl_2(array_hdl, int ndim1, int ndim2, - * int ndim1Old, int ndim2Old, const double val) - * - * Reallocates a two dimensional array of doubles. - * This routine will then copy the pertinent information from - * the old array to the new array. - * - * If both old dimensions are set to zero or less, then this routine - * will free the old memory before mallocing the new memory. This may - * be a benefit for extremely large mallocs. - * In all other cases, the new and the old malloced arrays will - * exist for a short time together. - * - * Input - * ------- - * array_hdl = Pointer to the global variable that - * holds the old and (eventually new) - * address of the array of doubles to be reallocated - * ndim1 = First dimension of the new array - * ndim2 = Second dimension of the new array - * ndim1Old = First dimension of the old array - * ndim2Old = Second dimension of the old array - * val = Default fill value. - **************************************************************************/ - { - if (ndim1 <= 0) ndim1 = 1; - if (ndim2 <= 0) ndim2 = 1; - ndim1Old = MaxI(ndim1Old, 0); - ndim2Old = MaxI(ndim2Old, 0); - /* - * One way to do it, if old information isn't needed. In this algorithm - * the arrays are never malloced at the same time. - */ - if ((*array_hdl == NULL) || (ndim1Old <= 0 && ndim2Old <= 0)) { - mdp_safe_free((void **) array_hdl); - *array_hdl = mdp_alloc_dbl_2(ndim1, ndim2, val); - if (*array_hdl == NULL) mdp_alloc_eh2("mdp_realloc_dbl_2"); - } - /* - * Other way to do when old information is available and needed - */ - else { - double **array_old = *array_hdl; - *array_hdl = (double **) mdp_array_alloc(2, ndim1, ndim2, sizeof(double)); - if (*array_hdl == NULL) { - mdp_alloc_eh2("mdp_realloc_dbl_2"); - } else { - /* - * Now, let's initialize the arrays - */ - int ndim1Min = MinI(ndim1, ndim1Old); - int ndim2Min = MinI(ndim2, ndim2Old); - double **array_new = *array_hdl; - /* - * When the second dimensions are equal, we can copy blocks - * using the very efficient bit moving kernels. - */ - if (ndim2 == ndim2Old) { - size_t sz = ndim1Min * ndim2 * sizeof(double); - (void) memcpy((void *) array_new[0], (const void *) array_old[0], sz); - } - /* - * If the second dimensions aren't equal, then we have to - * break up the bit operations even more - */ - else { - size_t sz = ndim2Min * sizeof(double); - size_t sz2 = (ndim2 - ndim2Min) * sizeof(double); - for (int i = 0; i < ndim1Min; i++) { - (void) memcpy((void *) array_new[i], - (const void *) array_old[i], sz); - if (ndim2 > ndim2Min && val != MDP_DBL_NOINIT) { - if (val == 0.0) { - (void) memset((void *) (array_new[i] + ndim2Min), 0, sz2); - } else { - double *dptr = array_new[i]; - for (int j = ndim2Min; j < ndim2; j++) dptr[j] = val; - } - } - } - } - /* - * finish up initializing the rest of the array - */ - if (ndim1 > ndim1Min && val != MDP_DBL_NOINIT) { - if (val == 0.0) { - size_t sz = (ndim1 - ndim1Min) * ndim2 * sizeof(double); - (void) memset((void *) array_new[ndim1Min], 0, sz); - } else { - double *dptr = array_new[ndim1Min]; - int num = (ndim1 - ndim1Min) * ndim2; - for (int i = 0; i < num; i++) { - dptr[i] = val; - } - } - } - /* - * Free the old array - */ - mdp_safe_free((void **) &array_old); - } - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - char **mdp_alloc_VecFixedStrings(int numStrings, int lenString) - - /************************************************************************* - * - * mdp_alloc_VecFixedStrings: - * - * Allocate and initialize a vector of fixed-length - * strings. Each string is initialized to the NULL string. - * - * Input - * ------- - * numStrings = Number of strings - * lenString = Length of each string including the trailing null - * character - * Return - * ------ - * This value is initialized to the correct address of the array. - * A NULL value in the position indicates an error. - **************************************************************************/ - { - int i; - char **array; - if (numStrings <= 0) numStrings = 1; - if (lenString <= 0) lenString = 1; - array = (char **) mdp_array_alloc(2, numStrings, lenString, sizeof(char)); - if (array != NULL) { - for (i = 0; i < numStrings; i++) { - array[i][0] = '\0'; - array[i][lenString - 1] = '\0'; - } - } else { - mdp_alloc_eh("mdp_alloc_VecFixedStrings", - sizeof(char) * numStrings * lenString + - numStrings * sizeof(void *)); - } - return array; - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_realloc_VecFixedStrings(char ***array_hdl, int numStrings, - int numOldStrings, int lenString) - - /************************************************************************* - * - * mdp_realloc_VecFixedStrings: - * - * Reallocate and initialize a vector of fixed-length - * strings. Each new string is initialized to the NULL string. - * old strings are copied. - * - * Input - * ------- - * ***array_hdl = The pointer to the char ** location holding - * the data to be reallocated. - * numStrings = Number of strings - * numOldStrings = Number of old strings - * lenString = Length of each string including the trailing null - * character - **************************************************************************/ - { - int i; - char **array, **ao; - if (numStrings <= 0) numStrings = 1; - if (numStrings == numOldStrings) return; - if (lenString <= 0) lenString = 1; - array = (char **) mdp_array_alloc(2, numStrings, lenString, sizeof(char)); - if (array != NULL) { - int len = MinI(numStrings, numOldStrings); - ao = *array_hdl; - if (ao) { - for (i = 0; i < len; i++) { - strncpy(array[i], ao[i], lenString); - } - } - if (numStrings > numOldStrings) { - for (i = numOldStrings; i < numStrings; i++) { - array[i][0] = '\0'; - array[i][lenString - 1] = '\0'; - } - } - mdp_safe_free((void **) array_hdl); - *array_hdl = array; - - } else { - mdp_alloc_eh("mdp_realloc_VecFixedStrings", - sizeof(char) * numStrings * lenString + - numStrings * sizeof(void *)); - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_safe_alloc_VecFixedStrings(char ***array_hdl, - int numStrings, int lenString) - - /************************************************************************* - * - * mdp_safe_alloc_VecFixedStrings - * - * Allocate and initialize an array of strings of fixed length - * - * Input - * ------- - * *array_hdl = Previous value of pointer. If non-NULL will try - * to free the memory at this address. - * numStrings = Number of strings - * lenString = Length of each string including the trailing null - * character - * Output - * ------ - * *array_hdl = This value is initialized to the correct address - * of the array. - * A NULL value in the position indicates an error. - **************************************************************************/ - { - if (array_hdl == NULL) { - mdp_alloc_eh("mdp_safe_alloc_VecFixedStrings: handle is NULL", - MDP_ALLOC_INTERFACE_ERROR); - return; - } - if (*array_hdl != NULL) mdp_safe_free((void **) array_hdl); - *array_hdl = mdp_alloc_VecFixedStrings(numStrings, lenString); - if (*array_hdl == NULL) - mdp_alloc_eh2("mdp_safe_alloc_VecFixedStrings"); - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void **mdp_alloc_ptr_1(int numPointers) - - /************************************************************************* - * - * mdp_alloc_ptr_1: - * - * Allocate and initialize a vector of pointers - * of type pointer to void. All pointers are initialized to the NULL - * value. - * - * Input - * ------- - * numPointers = Number of pointers - * Return - * ------ - * This value is initialized to the correct address of the vector. - * A NULL value in the position indicates an error. - **************************************************************************/ - { - int i; - void **array; - if (numPointers <= 0) numPointers = 1; - array = (void **) mdp_array_alloc(1, numPointers, sizeof(void *)); - if (array != NULL) { - for (i = 0; i < numPointers; i++) { - array[i] = NULL; - } - } else { - mdp_alloc_eh("mdp_alloc_ptr_1", - sizeof(void *) * numPointers); - } - return array; - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_safe_alloc_ptr_1(void ***array_hdl, int numPointers) - - /************************************************************************** - * - * mdp_safe_alloc_ptr_1: - * - * Allocate and initialize a vector of pointers - * of type pointer to void. All pointers are initialized to the NULL - * value. - * - * Input - * ------- - * *array_hdl = Previous value of pointer. If non-NULL will try - * to free the memory at this address. - * numPointers = Number of pointers - * Output - * ------ - * *array_hdl = This value is initialized to the correct address - * of the array. - * A NULL value in the position indicates an error. - **************************************************************************/ - { - if (array_hdl == NULL) { - mdp_alloc_eh("mdp_safe_alloc_ptr_1: handle is NULL", - MDP_ALLOC_INTERFACE_ERROR); - return; - } - if (*array_hdl != NULL) mdp_safe_free((void **) array_hdl); - *array_hdl = mdp_alloc_ptr_1(numPointers); - if (*array_hdl == NULL) mdp_alloc_eh2("mdp_safe_alloc_ptr_1"); - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_realloc_ptr_1(void ***array_hdl, int numLen, int numOldLen) - - /************************************************************************* - * - * mdp_realloc__ptr_1: - * - * Reallocate and initialize a vector of pointers - * Each new pointer is initialized to NULL. - * old Pointers are copied. - * - * Input - * ------- - * ***array_hdl = The pointer to the char ** location holding - * the data to be reallocated. - * numLen = Number of strings - * numOldLen = Number of old strings - **************************************************************************/ - { - if (array_hdl == NULL) { - mdp_alloc_eh("mdp_safe_alloc_ptr_1: handle is NULL", - MDP_ALLOC_INTERFACE_ERROR); - return; - } - if (numLen <= 0) numLen = 1; - if (numOldLen < 0) numOldLen = 0; - if (numLen == numOldLen) return; - size_t bytenum = sizeof(void *) * numLen; - void **array = (void **) smalloc(bytenum); - if (array != NULL) { - int len = MinI(numLen, numOldLen); - if (*array_hdl) { - void **ao = *array_hdl; - for (int i = 0; i < len; i++) array[i] = ao[i]; - } else { - numOldLen = 0; - } - if (numLen > numOldLen) { - bytenum = sizeof(void *) * (numLen - numOldLen); - (void) memset((void *) (array + numOldLen), 0, bytenum); - } - mdp_safe_free((void **) array_hdl); - *array_hdl = array; - } else { - mdp_alloc_eh("mdp_realloc_ptr_1", sizeof(void *) * numLen); - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_copy_VecFixedStrings(char ** const copyTo, - const char ** const copyFrom, - int numStrings, size_t maxLenString) - - /************************************************************************* - * - * mdp_copy_VecFixedStrings - * - * Copies an array of string vectors. - * - * Input - * ------- - * copyFrom = vector of C strings. It should be null terminated - * numStrings = number of strings - * maxLenString = maximum of the size of the string arrays, - * copyTo and copyFrom. This is used as the - * argument to strncpy() function. - * - * Output - * copyTo = vector of strings - * - **************************************************************************/ - { - if (maxLenString > 0) { - char *dTo; - const char *dFrom; - for (int i = 0; i < numStrings; i++) { - dTo = copyTo[i]; - dFrom = copyFrom[i]; - if (!dTo) { - mdp_alloc_eh2("mdp_copy_VecFixedStrings"); - } - if (!dFrom) { - mdp_alloc_eh2("mdp_copy_VecFixedStrings"); - } - (void) strncpy(dTo, dFrom, maxLenString); - dTo[maxLenString-1] = '\0'; - } - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - char *mdp_copy_string(const char * const copyFrom) - - /************************************************************************* - * - * mdp_copy_string: - * - * Allocates space for and copies a string - * - * Input - * ------- - * copyFrom = null terminated string. If NULL is supplied, then - * nothing is malloced and a NULL value is returned. - * Return - * ------ - * This value is initialized to the correct address of the array. - * A NULL value in the position either indicates an error, or - * that the original pointer to the string was NULL. - **************************************************************************/ - { - char *cptr; - if (copyFrom == NULL) return NULL; - cptr = (char *) mdp_array_alloc(1, strlen(copyFrom) + 1, sizeof(char)); - if (cptr != NULL) { - (void) strcpy(cptr, copyFrom); - } else { - mdp_alloc_eh("mdp_copy_string", - static_cast(sizeof(char) * (strlen(copyFrom) + 1))); - } - return cptr; - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_safe_copy_string(char **string_hdl, const char *copyFrom) - - /************************************************************************* - * - * mdp_safe_copy_string: - * - * Allocates space for and copies a string - * - * Input - * ------- - * *string_hdl = Previous value of pointer. If non-NULL will try - * to free the memory at this address. - * *copyFrom = String to be copied - * Output - * ------ - * *string_hdl = Pointer to the copied string - * A NULL value in the position indicates an error. - **************************************************************************/ - { - if (string_hdl == NULL) { - mdp_alloc_eh("mdp_safe_copy_string: string_hdl is NULL", - MDP_ALLOC_INTERFACE_ERROR); - return; - } - if (*string_hdl != NULL) mdp_safe_free((void **) string_hdl); - if (copyFrom == NULL) { - *string_hdl = NULL; - return; - } - *string_hdl = mdp_copy_string(copyFrom); - if (*string_hdl == NULL) { - mdp_alloc_eh2("mdp_safe_copy_string"); - } - return; - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_copy_dbl_1(double * const copyTo, - const double * const copyFrom, const int len) - - /************************************************************************** - * - * mdp_copy_dbl_1: - * - * Copies one contiguous double vector into another double vector - * - * Input - * ------------- - * *copyFrom = Vector of values to be copied - * len = Length of the vector - * - * Output - * ------------ - * *copyTo = Vector of values to receive the copy - * - **************************************************************************/ - { - if (len > 0) { - (void) memcpy((void *)copyTo, (const void *)copyFrom, len * sizeof(double)); - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_copy_dbl_2(double ** const copyTo, - const double ** const copyFrom, const int len1, - const int len2) - - /************************************************************************** - * - * mdp_copy_dbl_2: - * - * Copies one double array into another double array - * - * Input - * ------------- - * *copyFrom = Vector of values to be copied - * len1 = Length of the first array - * len2 = length of the second array - * - * Output - * ------------ - * *copyTo = array of values to receive the copy - * - **************************************************************************/ - { - if (len1 > 0 && len2 > 0) { - size_t bytelen = len1 * len2 * sizeof(double); - const double * const dFrom = copyFrom[0]; - double * const dTo = copyTo[0]; - (void) memcpy((void *)dTo, (const void *)dFrom, bytelen); - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_copy_int_1(int *const copyTo, - const int * const copyFrom, int len) - - /************************************************************************** - * - * mdp_copy_int_1: - * - * Copies one int vector into int double vector - * - * Input - * ------------- - * *copyFrom = Vector of values to be copied - * len = Length of the vector - * - * Output - * ------------ - * *copyTo = Vector of values to receive the copy - * - **************************************************************************/ - { - if (len > 0) { - size_t bytelen = len * sizeof(int); - (void) memcpy((void *)copyTo, (const void *)copyFrom, bytelen); - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_copy_int_2(int **const copyTo, - const int ** const copyFrom, const int len1, - const int len2) - - /************************************************************************** - * - * mdp_copy_int_2: - * - * Copies one 2D int array into another 2D int array - * - * Input - * ------------- - * *copyFrom = Vector of values to be copied - * len1 = Length of the first array - * len2 = Length of the second array - * - * Output - * ------------ - * *copyTo = array of values to receive the copy - * - **************************************************************************/ - { - if (len1 > 0 && len2 > 0) { - size_t bytelen = len1 * len2 * sizeof(int); - const int * const dFrom = copyFrom[0]; - int * const dTo = copyTo[0]; - (void) memcpy((void *)dTo, (const void *)dFrom, bytelen); - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_copy_ptr_1(void *const copyTo, - const void * const copyFrom, int len) - - /************************************************************************** - * - * mdp_copy_ptr_1: - * - * Copies one ptr vector into another ptr vector - * - * Input - * ------------- - * *copyFrom = Vector of ptr values to be copied - * len = Length of the vector - * - * Output - * ------------ - * *copyTo = Vector of values to receive the copy - * - **************************************************************************/ - { - if (len > 0) { - size_t bytelen = len * sizeof(void *); - (void) memcpy((void *)copyTo, (const void *)copyFrom, bytelen); - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void **mdp_dupl_ptr_1(const void * const copyFrom, int len) - - /************************************************************************** - * - * mdp_dupl_ptr_1: - * - * duplicates one ptr vector into another ptr vector - * - * Input - * ------------- - * *copyFrom = Vector of ptr values to be copied - * len = Length of the vector - * - * Output - * ------------ - * *copyTo = Vector of values to receive the copy - * - **************************************************************************/ - { - if (len > 0) { - void **array = mdp_alloc_ptr_1(len); - if (copyFrom) { - size_t bytelen = len * sizeof(void *); - (void) memcpy((void *)array, (const void *)copyFrom, bytelen); - } - return array; - } - return 0; - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_init_dbl_1(double * const v, const double value, const int len) - - /************************************************************************** - * - * mdp_init_dbl_1: - * - * Assigns a single value to a double vector - * - * Input - * ------------- - * v = Vector of values to be assigned - * value = value to assign with - * len = Length of the vector - * - **************************************************************************/ - { - if (len > 0) { - if (value == 0.0) { - (void) memset((void *)v, 0, len * sizeof(double)); - } else { - int m = len % 7; - if (m != 0) { - for (int i = 0; i < m; m++) { - v[i] = value; - } - if (len < 7) return; - } - for (int i = m; i < len; i += 7) { - v[i ] = value; - v[i+1] = value; - v[i+2] = value; - v[i+3] = value; - v[i+4] = value; - v[i+5] = value; - v[i+6] = value; - } - } - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_zero_dbl_1(double * const v, const int len) - - /************************************************************************** - * - * mdp_zero_dbl_1: - * - * Zeroes out a double vector (special form of mdp_allo_dbl_1()) - * - * Input - * ------------- - * v = Vector of values to be set to zero - * len = Length of the vector - **************************************************************************/ - { - if (len > 0) { - (void) memset((void *)v, 0, len * sizeof(double)); - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_init_dbl_2(double ** const v, double value, int len1, int len2) - - /************************************************************************** - * - * mdp_init_dbl_2: - * - * Assigns a single value to a double matrix. Contiguous data for the - * matrix is assumed. - * - * Input - * ------------- - * v = matrix of values to be assigned - * value = value to assign with - * len = Length of the vector - * - **************************************************************************/ - { - int len = len1 * len2; - if (len > 0 && len1 > 0 && v) { - double * const dstart = v[0]; - if (value == 0.0) { - size_t bytelen = len * sizeof(double); - (void) memset((void *)dstart, 0, bytelen); - } else { - int m = len % 7; - if (m != 0) { - for (int i = 0; i < m; m++) { - dstart[i] = value; - } - if (len < 7) return; - } - for (int i = m; i < len; i += 7) { - dstart[i ] = value; - dstart[i+1] = value; - dstart[i+2] = value; - dstart[i+3] = value; - dstart[i+4] = value; - dstart[i+5] = value; - dstart[i+6] = value; - } - } - } - } - /****************************************************************************/ - /****************************************************************************/ - /****************************************************************************/ - - void mdp_init_int_1(int * const v, const int value, const int len) - - /************************************************************************** - * - * mdp_init_int_1: - * - * Assigns a single value to an int vector - * - * Input - * ------------- - * v = Vector of values to be assigned - * value = value to assign with - * len = Length of the vector - * - **************************************************************************/ - { - if (len > 0) { - if (value == 0) { - size_t bytelen = len * sizeof(int); - (void) memset((void *)v, 0, bytelen); - } else { - int m = len % 7; - if (m != 0) { - for (int i = 0; i < m; m++) { - v[i] = value; - } - if (len < 7) return; - } - for (int i = m; i < len; i += 7) { - v[i ] = value; - v[i+1] = value; - v[i+2] = value; - v[i+3] = value; - v[i+4] = value; - v[i+5] = value; - v[i+6] = value; - } - } - } - } -} -/****************************************************************************/ -/* END of mdp_allo.cpp */ -/****************************************************************************/ diff --git a/Cantera/src/base/mdp_allo.h b/Cantera/src/base/mdp_allo.h deleted file mode 100644 index 027c268a4..000000000 --- a/Cantera/src/base/mdp_allo.h +++ /dev/null @@ -1,734 +0,0 @@ -/** - * @file mdp_allo.h - * Declarations for Multi Dimensional Pointer (mdp) malloc routines, which - * allow for dimensioning of arbitrarily dimensioned pointer arrays using - * one call. - */ -/* - * Copywrite 2004 Sandia Corporation. Under the terms of Contract - * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government - * retains certain rights in this software. - * See file License.txt for licensing information. - */ - -#ifndef MDP_ALLO_H -#define MDP_ALLO_H - -#include - -/* - * Include the header here in order to pick up size_t definition - */ -#include - -/** - * The mdp routines are extremely lightweight and fast fortran compatibile - * malloc routines for allocating multiple dimensioned arrays of doubles - * ints, char, and pointers using a single call. These routines don't - * use the std C+ lib routines. - * - * All calls are essentially wrappers around the routine mdp_alloc_array() - * which allocates multidimensioned arrays. The arrays contain room for - * the data and the pointer information that is used to access data - * in the object. - * - * One convention that is always used is that a pointer that is - * not malloced always has a value of zero. If the pointer is nonnull, - * then it may be freed, always (and vica-versa). - * - * Where possible, the low leve routines - * memcopy and memset are used to copy or zero memory. - * - * No array bounds checking is ever done within these routines. buyer beware. - * The bounds of arrays are not carried with the array object, ever. - * Bounds of arrays are input via the parameter list. Normally, - * for applications this is not a problem, since the application - * knows what the array dimensions are. - * - * There are several other general principles. - * - * If an allocation size is set to 0, then the actual - * allocation size is set to 1 within the program. - * Something is always allocated whenever a call is made - * to an mdp routine. - * - * All checks for allocations are always checked for success. - * If a failure is found, thene mdp_alloc_eh() is called - * for disposition of the error condition. - * - * Error handling behavior is set by the MDP_ALLO_errorOption external - * int. The default error behavior is to print an error message to stderr, and - * then throw an exception that inherited from std::exception. Usually - * std::bad_alloc() is thrown whenever there is a problem. - * - */ -namespace mdp { - - /** - * If we have array_alloc() from another Sandia program, we will not use - * the one from this mdp_array_alloc. Instead we will redefine the names - */ -#ifdef HAVE_ARRAY_ALLOC -# define mdp_array_alloc array_alloc -# define mdp_safe_free safe_free -#endif - - /*! - * MDP_INT_NOINIT is a poor man's way of specifying whether a value should be - * initialized. These are seldom used numbers which can be used in place - * of real ints and dbls to indicate that initialization shouldn't take - * place. - */ -#define MDP_INT_NOINIT -68361 - - /*! - * MDP_DBL_NOINIT is a poor man's way of specifying whether a value should be - * initialized. These are seldom used numbers which can be used in place - * of real ints and dbls to indicate that initialization shouldn't take - * place. - */ -#define MDP_DBL_NOINIT -1.241E11 - - /*! - * Error Handling - * 7 print and exit - * 6 exit - * 5 print and create a divide by zero for stack trace analysis. - * 4 create a divide by zero for stack analysis trace - * 3 print a message and throw the std::bad_alloc() exception. - * 2 throw the std::bad_alloc() exception and be quiet. - * 1 print a message and return from package with the NULL pointer - * 0 Keep completely silent about the matter and return with - * a null pointer. - * - * -> Right now, the only way to change this option is to right here - * - * The default is to set it to 3. - */ - extern int MDP_ALLO_errorOption; - - - /****************************************************************************/ - /* - * Externals that should be set by the calling program. - * These are only used for debugging purposes. - */ -#ifdef MDP_MPDEBUGIO - extern int MDP_MP_Nprocs; - extern int MDP_MP_myproc; -#endif - - /****************************************************************************/ - /* - * MDP_SAFE_DELETE() - * This useful define is great for delete single instances of - * mallocing using new. - */ -#define MDP_SAFE_DELETE(a) if (a) { delete (a); a = 0; } - - - /****************************************************************************/ - -#define mdp_alloc_struct(x, num) (x *) mdp_array_alloc(1, (num), sizeof(x)) - - - /* function declarations for dynamic array allocation */ - - //! allocates multidimensional pointer arrays of arbitrary length - //! via a single malloc call - /*! - * The first dimension is the number of dimensions in the allocation - */ - extern double *mdp_array_alloc(int numdim, ...); - - //! Free a vector and set its value to 0 - /*! - * This function carries out the following operation - * @code - * free(*hndVec); - * *hndVec = 0; - * @endcode - * - * @param hndVec This is the address of the pointer, expressed as - * a void ** - * - * Note, a key idea behind the mdp suite of routines is that - * a pointer that can be malloced is either malloced or its value - * is 0. This routine enforces this convention. - */ - extern void mdp_safe_free(void ** hndVec); - - //! Allocate a vector of integers - /*! - * The vector is initialized, unless the default int value is set - * to MDP_INT_NOINIT - * - * @param len Length of the vector - * @param defval Default value for the int, defaults to MDP_INT_NOINIT - * - * @return returns a pointer to the vector - */ - extern int *mdp_alloc_int_1(int len, const int defval = MDP_INT_NOINIT); - - //! Allocate a vector of integers, potentially freeing memory first - /*! - * The vector is initialized, unless the default int value is set - * to MDP_INT_NOINIT - * - * Input - * -------------- - * @param array_hdl Previous value of pointer. If non-NULL will try - * to free the memory at this address before doing - * a new alloc - * @param len Length of the vector - * @param defval Default value for the int, defaults to MDP_INT_NOINIT - * - * Output - * --------- - * @return *array_hdl = This value is initialized to the correct address - * of the array. - * A NULL value in the position indicates an error. - */ - extern void mdp_safe_alloc_int_1(int **array_hdl, int len, - const int defval = MDP_INT_NOINIT); - - //! Reallocates a one dimensional array of ints, copying old - //! information to the new array - /*! - * Reallocates a one dimensional array of ints. - * This routine always allocates space for at least one int. - * Calls the smalloc() routine to ensure that all malloc - * calls go through one location. This routine will then copy - * the pertinent information from the old array to the - * new array. - * - * NewArray[0:old_len-1] = OldArray[0:old_len-1]; - * NewArray[old_len:new_len-1] = defval; - * - * Input - * -------------- - * @param array_hdl Previous value of pointer. If non-NULL will try - * to free the memory at this address before doing - * a new alloc - * @param new_len New Length of the vector - * @param old_len New Length of the vector - * @param defval Default value for the int, defaults to MDP_INT_NOINIT - * - * Output - * --------- - * @return *array_hdl = This value is initialized to the correct address - * of the array. - * A NULL value in the position indicates an error. - */ - extern void mdp_realloc_int_1(int **array_hdl, int new_len, int old_len, - const int defval = MDP_INT_NOINIT); - - - //! Allocate a 2D matrix of integers - /*! - * The matrix is initialized, unless the default int value is set - * to MDP_INT_NOINIT, which is the default. - * - * matrix[len1][len2] - * - * All int data entries are contiguous. Therefore, it may - * be used as input into BLAS matrix function calls. - * This can be considered to be in fortran order format with - * len2 as the number of rows, and len1 as the number of columns. - * - * matrix[jcol] refers to the jcol column of the matrix. - * Therefore, matrix[0] is a pointer to the beginning of the - * data portion of the structure. - * The structure will have len1 pointers at the beginning - * that holds pointers into the top of the columns of the - * contiguous data. - * - * The entire structure may be deallocated via one free call. - * - * @param len1 Outer Length of the vector - * @param len2 Inner length of the matrix - * @param defval Default value for the int, defaults to MDP_INT_NOINIT - * - * @return returns a pointer to the matrix - */ - extern int **mdp_alloc_int_2(int len1, int len2, - const int defval = MDP_INT_NOINIT); - - - //! Allocate and initialize a one dimensional array of doubles. - /*! - * As per the convention in mdp, this routine always initializes - * at least one slot. - * - * @param nvalues Length of the array. If this number is - * less than one, it is set to one. Therefore, - * This routine always initializes at least - * one double. - * @param val intialization value. Set it to the - * constant MDP_DBL_NOINIT if you don't - * want any initialization. memset() is - * used for zero initialization for fast - * execution speed. - * - * @return Pointer to the intialized array of doubles - * Failures are indicated by returning the NULL pointer. - */ - extern double *mdp_alloc_dbl_1(int nvalues, const double val=MDP_DBL_NOINIT); - - //! Allocates and/or initialises a one dimensional array of doubles. - /*! - * This routine will free any old memory that was located at that - * position, before it will allocate a new vector. - * - * @param hndVec Previous value of pointer. If non-NULL will try - * to free the memory at this address. On output, - * this value is initialized to the correct address - * of the array. A NULL value in the position - * indicates an error. - * @param nvalues Length of the array - * @param val initialization value - * - */ - extern void mdp_safe_alloc_dbl_1(double **hndVec, int nvalues, - const double val=MDP_DBL_NOINIT); - - //! Reallocate a vector of doubles possibly retaining a subset of values - /*! - * Reallocates the array and sets: - * - * (*hndVec)[0:oldLen-1] = oldVec[0:oldLen-1] - * (*hndVec)[oldLen:newLen] = defVal - * - * Input - * ------ - * @param newLen New Length of the vector - * - * Output - * ------- - * @param oldLen Old Length of the vector - * - */ - extern void mdp_realloc_dbl_1(double ** hndVec, int newLen, int oldLen, - const double defVal=MDP_DBL_NOINIT); - - //! Allocate and initialize a two dimensional array of doubles. - /*! - * Allocate a two dimensional array of doubles. The array is in - * fortran order and can be accessed via the following form: - * - * dblArray[ndim1][ndim2] - * - * Note, ndim2 is the inner dimension. i.e., the array is - * in column ordering. - * - * Input - * ------- - * @param ndim1 Length of the first dimension of the array - * @param ndim2 Length of the second dimension of the array - * @param val Intialization value - * - * - * @return Pointer to the initialized array of doubles. - * Failures are indicated by returning the NULL pointer. - */ - extern double **mdp_alloc_dbl_2(int ndim1, int ndim2, const double val); - - //! Allocate and initialize a two dimensional array of doubles. - /*! - * Allocate a two dimensional array of doubles. The array is in - * fortran order and can be accessed via the following form: - * - * (*arrayHndl)[ndim1][ndim2] - * - * Note, ndim2 is the inner dimension. i.e., the array is - * in column ordering. - * - * Input - * ------- - * @param ndim1 Length of the first dimension of the array - * @param ndim2 Length of the second dimension of the array - * @param val Intialization value - * @param arrayHndl Handle to the array. If nonnull, the array - * is first freed. Failures are indicated - * by returning the NULL pointer. - */ - extern void mdp_safe_alloc_dbl_2(double ***arrayHndl, int ndim1, int ndim2, - const double val = MDP_DBL_NOINIT); - - //! Reallocates a two dimensional array of doubles to a new set of - //! dimensions, copying the old results into the new array. - /*! - * This routine will then copy the pertinent information from - * the old array to the new array. - * - * If both old dimensions are set to zero or less, then this routine - * will free the old memory before mallocing the new memory. This may - * be a benefit for extremely large mallocs. - * In all other cases, the new and the old malloced arrays will - * exist for a short time together. - * - * Input - * ------- - * @param hndArray = Pointer to the global variable that - * holds the old and (eventually new) - * address of the array of doubles to be reallocated - * @param ndim1 = First dimension of the new array - * @param ndim2 = Second dimension of the new array - * @param ndim1Old = First dimension of the old array - * @param ndim2Old = Second dimension of the old array - * @param defVal = Default fill value. - * - * Output - * -------------- - * The resulting vector looks like this: - * - * (*hndArray)[0:ndim1Old-1][0:ndim2Old-1] = oldVec[0:ndim1Old-1][0:ndim2Old-1] - * (*hndArray)[ndim1Old:ndim1][ndim2Old:ndim2] = defVal - * - */ - extern void mdp_realloc_dbl_2(double *** hndArray, int ndim1, int ndim2, - int ndim1Old, int ndim2Old, - const double defVal=MDP_DBL_NOINIT); - - //! Allocate and initialize a one dimensional array of characters. - /*! - * The array is always initialized. - * - * Input - * ------- - * @param nvalues Length of the array - * @param val intialization value. defaults to the NULL char - * - * @return Pointer to the intialized character array - * Failures are indicated by returning the NULL pointer. - */ - extern char *mdp_alloc_char_1(int nvalues, const char val = '\0'); - - //! Allocate and initialize a one dimensional array of characters, - //! deallocating the space before hand. - /*! - * This routine will free any old memory that was located at that - * position, before it will allocate a new vector. - * - * The array is always initialized. - * - * @param arrayHnd Pointer to the global variable that - * holds the old and (eventually new) - * address of the array of char to be reallocated - * @param nvalues Length of the array - * @param val intialization value. defaults to the NULL char - * - * @return Pointer to the intialized character array - * Failures are indicated by returning the NULL pointer. - */ - extern void mdp_safe_alloc_char_1(char **arrayHnd, int nvalues, - const char val = '\0'); - - //! Allocate and initialize a vector of fixed-length - //! strings. Each string is initialized to the NULL string. - /*! - * @param numStrings Number of strings - * @param lenString Length of each string including the trailing null - * character - * - * @return Value is initialized to the correct address - * of the new array on exit. - * A NULL value in the position indicates an error. - */ - extern char **mdp_alloc_VecFixedStrings(int numStrings, int lenString); - - //! Allocate and initialize an array of strings of fixed length - /*! - * @param numStrings Number of strings - * @param lenString Length of each string including the trailing null - * character - * - * @param array_hdl Value is initialized to the correct address - * of the new array on exit. - * A NULL value in the position indicates an error. - * If non-NULL on entry the routine will first - * free the memory at the address. - */ - extern void mdp_safe_alloc_VecFixedStrings(char ***arrayHnd, int numStrings, - int lenString); - - //! Reallocate and initialize a vector of fixed-length strings. - /*! - * Each new string is initialized to the NULL string. - * Old strings are copied. - * - * @param array_hdl The pointer to the char ** location holding - * the data to be reallocated. - * @param numStrings Number of strings - * @param numOldStrings Number of old strings - * @param lenString Length of each string including the trailing null - * character - */ - extern void mdp_realloc_VecFixedStrings(char ***array_hdl, int numStrings, - int numOldStrings, int lenString); - - //! Allocate and initialize a vector of pointers of type pointer to void. - /*! - * All pointers are initialized to the NULL value. - * - * @param numPointers Number of pointers - * - * @return This value is initialized to the correct address of the vector. - * A NULL value in the position indicates an error. - */ - extern void **mdp_alloc_ptr_1(int numPointers); - - //! Allocate and initialize a vector of pointers of type pointer to void. - /*! - * All pointers are initialized to the NULL value. - * - * @param numPointers Number of pointers - * @param array_hdl This value is initialized to the correct address - * of the array. - * A NULL value in the position indicates an error. - * Previous value of pointer. If non-NULL will try - * to free the memory at this address. - */ - extern void mdp_safe_alloc_ptr_1(void ***array_hnd, int numPointers); - - //! Reallocate and initialize a vector of pointers - /*! - * All old pointers are copied - * Each new pointer not associated with an old pointer is - * initialized to NULL. - * - * @param array_hdl The pointer to the char ** location holding - * the data to be reallocated. - * @param numLen Number of new pointers - * @param numOldLen Number of old pointers - */ - extern void mdp_realloc_ptr_1(void ***array_hdl, int numLen, int numOldLen); - - //! Copies one ptr vector into another ptr vector - /*! - * - * @param copyFrom Vector of ptr values to be copied - * @param len Length of the vector - * - * @param copyTo Vector of values to receive the copy - */ - extern void mdp_copy_ptr_1(void *const copyTo, - const void * const copyFrom, const int len); - - //! Duplicates one ptr vector into another ptr vector - /*! - * Mallocs a copy of one vector of pointers and returns the pointer - * to the copy. - * - * Input - * ------------- - * @param *copyFrom Vector of ptr values to be copied - * @param len Length of the vector - * - * Output - * ------------ - * @return Vector of values to receive the copy - */ - extern void **mdp_dupl_ptr_1(const void * const copyFrom, int len); - - //! Copies an array of string vectors from one char ** vector to - //! another - /*! - * The space must have already been allocated within copyFrom and copyTo - * arrays. Overwrites are prevented by the proper application of the - * variable maxLenString. Strings are forced to be null-terminated - * Therefore copyTo[maxLenString-1] = '/0' - * - * Input - * ------- - * @param copyFrom vector of C strings. It should be null terminated - * @param numStrings number of strings - * @param maxLenString maximum of the size of the string arrays, - * copyTo and copyFrom. This is used as the - * argument to strncpy() function. - * - * Output - * ------ - * @param copyTo vector of strings - * - */ - extern void mdp_copy_VecFixedStrings(char ** const copyTo, - const char ** const copyFrom, - int numStrings, size_t maxLenString); - - //! Allocates space for and copies a string - /*! - * - * @param copyFrom null terminated string. If NULL is supplied, then - * nothing is malloced and a NULL value is returned. - * - * @return This value is initialized to the correct address of the array. - * A NULL value in the position either indicates an error, or - * that the original pointer to the string was NULL. - */ - extern char *mdp_copy_string(const char * const copyFrom); - - //! Allocates space for and copies a string - /*! - * @param stringHnd Previous value of pointer. If non-NULL will try - * to free the memory at this address. - * - * @param copyFrom null terminated string. If NULL is supplied, then - * nothing is malloced and a NULL value is returned. - * - * @return This value is initialized to the correct address of the array. - * A NULL value in the position either indicates an error, or - * that the original pointer to the string was NULL. - */ - extern void mdp_safe_copy_string(char **stringHnd, const char * const copyFrom); - - //! Copy a double vector to a double vector - /*! - * copyTo[len] = copyFrom[len] - * - * Input - * ------- - * @param copyFrom Vector to copy ( length >= len) - * @param len Length of the copy - * - * Output - * ------- - * @param copyTo Vector to receive the copy ( length >= len) - */ - extern void mdp_copy_dbl_1(double * const copyTo, - const double * const copyFrom, - const int len); - - //! Copy a double array to a double array - /*! - * This routine carries out a straight copy on the effective 1D description - * of each of the arrays. It copies - * the first len1*len2 doubless storred within copyFrom into the - * the first len1*len2 double slots in copyTo. It does not account - * for the actual dimensions of the two arrays. - * - * Input - * -------- - * @param copyFrom Vector to copy ( length >= len1 * len2) - * @param len1 Length of the first dimension - * @param len2 Length of the second dimension - * - * Output - * ---------- - * @param copyTo Array to receive the copy ( length >= len1 * len2) - */ - extern void mdp_copy_dbl_2(double ** const copyTo, const double ** const copyFrom, - const int len1, const int len2); - - //! Copies one int vector into one int vector - /*! - * Input - * ------------- - * @param copyFrom Vector of values to be copied - * @param len Length of the vector - * - * Output - * ------------ - * *copyTo = Vector of values to receive the copy - */ - extern void mdp_copy_int_1(int * const copyTo, const int * const copyFrom, - const int len); - - //! Copies one 2D int array into another 2D int array - /*! - * This routine carries out a straight copy. Actually it copies - * the first len1*len2 ints storred within copyFrom into the - * the first len1*len2 int slots in copyTo. It does not account - * for the actual dimensions of the two arrays. - * - * @param copyFrom Vector of values to be copied - * @param len1 Length of the first array - * @param len2 Length of the second array - * Output - * ------------ - * @param copyTo Vector of values to receive the copy - */ - extern void mdp_copy_int_2(int ** const copyTo, const int ** const copyFrom, - const int len1, const int len2); - - //! Assigns a single value to a double vector - /*! - * @param v Vector of values to be assigned - * @param value Value to assign with - * @param len Length of the vector - */ - extern void mdp_init_dbl_1(double * const v, const double value, const int len); - - //! Zeroes a double vector - /*! - * @param v = Vector of values to be assigned - * @param len = Length of the vector - */ - extern void mdp_zero_dbl_1(double * const v , const int len); - - //! Assigns a single value to a double matrix. Contiguous data for the - //! matrix is assumed. - /*! - * Input - * ------------- - * @param v matrix of values to be assigned - * @param value value to assign with - * @param len1 Length one of the vector - * @param len2 length two of the vector - */ - extern void mdp_init_dbl_2(double ** const v, const double value, - const int len1, const int len2); - - //! Assigns a single value to an int vector - /*! - * @param v Vector of values to be assigned - * @param value Value to assign with - * @param len Length of the vector - */ - extern void mdp_init_int_1(int * const v, const int value, const int len); - - - /* - * Utility routines to check that a number is finite - */ - - //! Utility routine to check to see that a number is neither zero - //! nor indefinite. - /*! - * This check can be used before using the number in a denominator. - * - * @param tmp number to be checked - */ - extern void checkZeroFinite(const double tmp); - - //! Utility routine to check to see that a number is finite. - /*! - * @param tmp number to be checked - */ - extern void checkFinite(const double tmp); - - //! Utility routine to link checkFinte() to fortran program - /*! - * This routine is accessible from fortran, usually - * - * @param tmp Pointer to the number to check - * - * @todo link it into the usual way Cantera handles Fortran calls - */ - extern "C" void checkfinite_(double *tmp); - - //! utility routine to check that a double stays bounded - /*! - * This routine checks to see if a number stays bounded. The absolute - * value of the number is required to stay below the trigger. - * - * @param tmp Number to be checked - * @param trigger bounds on the number. Defaults to 1.0E20 - */ - extern void checkMagnitude(const double tmp, const double trigger = 1.0E20); - -} /* end of mdp namespace */ -/****************************************************************************/ -#endif -/****************************************************************************/ - diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp deleted file mode 100644 index c75ffbe3c..000000000 --- a/Cantera/src/numerics/NonlinearSolver.cpp +++ /dev/null @@ -1,963 +0,0 @@ -/** - * - * @file NonlinearSolver.cpp - * - * Damped Newton solver for 0D and 1D problems - */ -/* - * Copywrite 2004 Sandia Corporation. Under the terms of Contract - * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government - * retains certain rights in this software. - * See file License.txt for licensing information. - */ - -#include - -#include "SquareMatrix.h" -#include "NonlinearSolver.h" - -#include "clockWC.h" -#include "vec_functions.h" -#include - -#include "mdp_allo.h" -#include - -extern void print_line(const char *, int); - -#include -#include -#include - - -#ifndef MAX -#define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) -#define MIN(x,y) (( (x) < (y) ) ? (x) : (y)) -#endif - -using namespace std; - -namespace Cantera { - - - - //----------------------------------------------------------- - // Constants - //----------------------------------------------------------- - - const double DampFactor = 4; - const int NDAMP = 7; - - //----------------------------------------------------------- - // Static Functions - //----------------------------------------------------------- - - static void print_line(const char *str, int n) { - for (int i = 0; i < n; i++) { - printf("%s", str); - } - printf("\n"); - } - - // Default constructor - /* - * @param func Residual and jacobian evaluator function object - */ - NonlinearSolver::NonlinearSolver(ResidJacEval *func) : - m_func(func), - neq_(0), - delta_t_n(-1.0), - m_nfe(0), - m_colScaling(0), - m_rowScaling(0), - m_numTotalLinearSolves(0), - m_numTotalNewtIts(0), - m_min_newt_its(0), - filterNewstep(0), - time_n(0.0), - m_matrixConditioning(0), - m_order(1), - rtol_(1.0E-3), - atolBase_(1.0E-10) - { - neq_ = m_func->nEquations(); - - m_ewt.resize(neq_, rtol_); - m_y_n.resize(neq_, 0.0); - m_y_nm1.resize(neq_, 0.0); - m_colScales.resize(neq_, 1.0); - m_rowScales.resize(neq_, 1.0); - m_resid.resize(neq_, 0.0); - atolk_.resize(neq_, atolBase_); - doublereal hb = std::numeric_limits::max(); - m_y_high_bounds.resize(neq_, hb); - m_y_low_bounds.resize(neq_, -hb); - - for (int i = 0; i < neq_; i++) { - atolk_[i] = atolBase_; - m_ewt[i] = atolk_[i]; - } - } - - NonlinearSolver::NonlinearSolver(const NonlinearSolver &right) { - *this =operator=(right); - } - - - NonlinearSolver::~NonlinearSolver() { - } - - NonlinearSolver& NonlinearSolver::operator=(const NonlinearSolver &right) { - if (this == &right) { - return *this; - } - // rely on the ResidJacEval duplMyselfAsresidJacEval() function to - // create a deep copy - m_func = right.m_func->duplMyselfAsResidJacEval(); - - neq_ = right.neq_; - m_ewt = right.m_ewt; - m_y_n = right.m_y_n; - m_y_nm1 = right.m_y_nm1; - m_colScales = right.m_colScales; - m_rowScales = right.m_rowScales; - m_resid = right.m_resid; - m_y_high_bounds = right.m_y_high_bounds; - m_y_low_bounds = right.m_y_low_bounds; - delta_t_n = right.delta_t_n; - m_nfe = right.m_nfe; - m_colScaling = right.m_colScaling; - m_rowScaling = right.m_rowScaling; - m_numTotalLinearSolves = right.m_numTotalLinearSolves; - m_numTotalNewtIts = right.m_numTotalNewtIts; - m_min_newt_its = right.m_min_newt_its; - filterNewstep = right.filterNewstep; - time_n = right.time_n; - m_matrixConditioning = right.m_matrixConditioning; - m_order = right.m_order; - rtol_ = right.rtol_; - atolBase_ = right.atolBase_; - atolk_ = right.atolk_; - - return *this; - } - - // Create solution weights for convergence criteria - /* - * We create soln weights from the following formula - * - * wt[i] = rtol * abs(y[i]) + atol[i] - * - * The program always assumes that atol is specific - * to the solution component - * - * param y vector of the current solution values - */ - void NonlinearSolver::createSolnWeights(const double * const y) { - for (int i = 0; i < neq_; i++) { - m_ewt[i] = rtol_ * fabs(y[i]) + atolk_[i]; - } - } - - // set bounds constraints for all variables in the problem - /* - * - * @param y_low_bounds Vector of lower bounds - * @param y_high_bounds Vector of high bounds - */ - void NonlinearSolver::setBoundsConstraints(const double * const y_low_bounds, - const double * const y_high_bounds) { - for (int i = 0; i < neq_; i++) { - m_y_low_bounds[i] = y_low_bounds[i]; - m_y_high_bounds[i] = y_high_bounds[i]; - } - } - - /** - * L2 Norm of a delta in the solution - * - * The second argument has a default of false. However, - * if true, then a table of the largest values is printed - * out to standard output. - */ - double NonlinearSolver::solnErrorNorm(const double * const delta_y, - bool printLargest) - { - int i; - double sum_norm = 0.0, error; - for (i = 0; i < neq_; i++) { - error = delta_y[i] / m_ewt[i]; - sum_norm += (error * error); - } - sum_norm = sqrt(sum_norm / neq_); - if (printLargest) { - const int num_entries = 8; - double dmax1, normContrib; - int j; - int *imax = mdp::mdp_alloc_int_1(num_entries, -1); - printf("\t\tPrintout of Largest Contributors to norm " - "of value (%g)\n", sum_norm); - printf("\t\t I ysoln deltaY weightY " - "Error_Norm**2\n"); - printf("\t\t "); print_line("-", 80); - for (int jnum = 0; jnum < num_entries; jnum++) { - dmax1 = -1.0; - for (i = 0; i < neq_; i++) { - bool used = false; - for (j = 0; j < jnum; j++) { - if (imax[j] == i) used = true; - } - if (!used) { - error = delta_y[i] / m_ewt[i]; - normContrib = sqrt(error * error); - if (normContrib > dmax1) { - imax[jnum] = i; - dmax1 = normContrib; - } - } - } - i = imax[jnum]; - if (i >= 0) { - printf("\t\t %4d %12.4e %12.4e %12.4e %12.4e\n", - i, m_y_n[i], delta_y[i], m_ewt[i], dmax1); - } - } - printf("\t\t "); print_line("-", 80); - mdp::mdp_safe_free((void **) &imax); - } - return sum_norm; - } - - /** - * L2 Norm of the residual - * - * The second argument has a default of false. However, - * if true, then a table of the largest values is printed - * out to standard output. - */ - double NonlinearSolver::residErrorNorm(const double * const resid, - bool printLargest) - { - int i; - double sum_norm = 0.0, error; - for (i = 0; i < neq_; i++) { - error = resid[i] / m_rowScales[i]; - sum_norm += (error * error); - } - sum_norm = sqrt(sum_norm / neq_); - if (printLargest) { - const int num_entries = 8; - double dmax1, normContrib; - int j; - int *imax = mdp::mdp_alloc_int_1(num_entries, -1); - printf("\t\tPrintout of Largest Contributors to norm " - "of Residual (%g)\n", sum_norm); - printf("\t\t I resid rowScale weightN " - "Error_Norm**2\n"); - printf("\t\t "); print_line("-", 80); - for (int jnum = 0; jnum < num_entries; jnum++) { - dmax1 = -1.0; - for (i = 0; i < neq_; i++) { - bool used = false; - for (j = 0; j < jnum; j++) { - if (imax[j] == i) used = true; - } - if (!used) { - error = resid[i] / m_rowScales[i]; - normContrib = sqrt(error * error); - if (normContrib > dmax1) { - imax[jnum] = i; - dmax1 = normContrib; - } - } - } - i = imax[jnum]; - if (i >= 0) { - printf("\t\t %4d %12.4e %12.4e %12.4e \n", - i, resid[i], m_rowScales[i], normContrib); - } - } - printf("\t\t "); print_line("-", 80); - mdp::mdp_safe_free((void **) &imax); - } - return sum_norm; - } - - - /** - * setColumnScales(): - * - * Set the column scaling vector at the current time - */ - void NonlinearSolver::setColumnScales() { - m_func->calcSolnScales(time_n, DATA_PTR(m_y_n), DATA_PTR(m_y_nm1), - DATA_PTR(m_colScales)); - } - - - void NonlinearSolver::doResidualCalc(const double time_curr, const int typeCalc, - const double * const y_curr, - const double * const ydot_curr, double* const residual, - int loglevel) - { - - - // Calculate the current residual - // Put the current residual into the vector, delta_y[] - // We need to pull this out of this function and carry it in. - m_func->evalResidNJ(time_curr, delta_t_n, y_curr, ydot_curr, residual); - m_nfe++; - } - - - // Compute the undamped Newton step - /* - * Compute the undamped Newton step. The residual function is - * evaluated at the current time, t_n, at the current values of the - * solution vector, m_y_n, and the solution time derivative, m_ydot_n. - * The Jacobian is not recomputed. - * - * A factored jacobian is reused, if available. If a factored jacobian - * is not available, then the jacobian is factored. Before factoring, - * the jacobian is row and column-scaled. Column scaling is not - * recomputed. The row scales are recomputed here, after column - * scaling has been implemented. - */ - void NonlinearSolver::doNewtonSolve(const double time_curr, const double * const y_curr, - const double * const ydot_curr, double* const delta_y, - SquareMatrix& jac, int loglevel) - { - int irow, jcol; - - - //! multiply the residual by -1 - for (int n = 0; n < neq_; n++) { - delta_y[n] = -delta_y[n]; - } - - - /* - * Column scaling -> We scale the columns of the Jacobian - * by the nominal important change in the solution vector - */ - if (m_colScaling) { - if (!jac.m_factored) { - /* - * Go get new scales -> Took this out of this inner loop. - * Needs to be done at a larger scale. - */ - // setColumnScales(); - - /* - * Scale the new Jacobian - */ - double *jptr = &(*(jac.begin())); - for (jcol = 0; jcol < neq_; jcol++) { - for (irow = 0; irow < neq_; irow++) { - *jptr *= m_colScales[jcol]; - jptr++; - } - } - } - } - - // if (m_matrixConditioning) { - // if (jac.m_factored) { - // m_func->matrixConditioning(0, neq_, delta_y); - // } else { - //double *jptr = &(*(jac.begin())); - // m_func->matrixConditioning(jptr, neq_, delta_y); - // } - //} - - /* - * row sum scaling -> Note, this is an unequivical success - * at keeping the small numbers well balanced and - * nonnegative. - */ - if (m_rowScaling) { - if (! jac.m_factored) { - /* - * Ok, this is ugly. jac.begin() returns an vector iterator - * to the first data location. - * Then &(*()) reverts it to a double *. - */ - double *jptr = &(*(jac.begin())); - for (irow = 0; irow < neq_; irow++) { - m_rowScales[irow] = 0.0; - } - for (jcol = 0; jcol < neq_; jcol++) { - for (irow = 0; irow < neq_; irow++) { - m_rowScales[irow] += fabs(*jptr); - jptr++; - } - } - - jptr = &(*(jac.begin())); - for (jcol = 0; jcol < neq_; jcol++) { - for (irow = 0; irow < neq_; irow++) { - *jptr /= m_rowScales[irow]; - jptr++; - } - } - } - for (irow = 0; irow < neq_; irow++) { - delta_y[irow] /= m_rowScales[irow]; - } - } - - - /* - * Solve the system -> This also involves inverting the - * matrix - */ - (void) jac.solve(delta_y); - - - /* - * reverse the column scaling if there was any. - */ - if (m_colScaling) { - for (irow = 0; irow < neq_; irow++) { - delta_y[irow] *= m_colScales[irow]; - } - } - -#ifdef DEBUG_JAC - if (printJacContributions) { - for (int iNum = 0; iNum < numRows; iNum++) { - if (iNum > 0) focusRow++; - double dsum = 0.0; - vector_fp& Jdata = jacBack.data(); - double dRow = Jdata[neq_ * focusRow + focusRow]; - printf("\n Details on delta_Y for row %d \n", focusRow); - printf(" Value before = %15.5e, delta = %15.5e," - "value after = %15.5e\n", y_curr[focusRow], - delta_y[focusRow], - y_curr[focusRow] + delta_y[focusRow]); - if (!freshJac) { - printf(" Old Jacobian\n"); - } - printf(" col delta_y aij " - "contrib \n"); - printf("--------------------------------------------------" - "---------------------------------------------\n"); - printf(" Res(%d) %15.5e %15.5e %15.5e (Res = %g)\n", - focusRow, delta_y[focusRow], - dRow, RRow[iNum] / dRow, RRow[iNum]); - dsum += RRow[iNum] / dRow; - for (int ii = 0; ii < neq_; ii++) { - if (ii != focusRow) { - double aij = Jdata[neq_ * ii + focusRow]; - double contrib = aij * delta_y[ii] * (-1.0) / dRow; - dsum += contrib; - if (fabs(contrib) > Pcutoff) { - printf("%6d %15.5e %15.5e %15.5e\n", ii, - delta_y[ii] , aij, contrib); - } - } - } - printf("--------------------------------------------------" - "---------------------------------------------\n"); - printf(" %15.5e %15.5e\n", - delta_y[focusRow], dsum); - } - } - -#endif - - m_numTotalLinearSolves++; - } - - /************************************************************************** - * - * boundStep(): - * - * Return the factor by which the undamped Newton step 'step0' - * must be multiplied in order to keep all solution components in - * all domains between their specified lower and upper bounds. - * Other bounds may be applied here as well. - * - * Currently the bounds are hard coded into this routine: - * - * Minimum value for all variables: - 0.01 * m_ewt[i] - * Maximum value = none. - * - * Thus, this means that all solution components are expected - * to be numerical greater than zero in the limit of time step - * truncation errors going to zero. - * - * Delta bounds: The idea behind these is that the Jacobian - * couldn't possibly be representative if the - * variable is changed by a lot. (true for - * nonlinear systems, false for linear systems) - * Maximum increase in variable in any one newton iteration: - * factor of 2 - * Maximum decrease in variable in any one newton iteration: - * factor of 5 - */ - double NonlinearSolver::boundStep(const double* const y, - const double* const step0, const int loglevel) { - int i, i_lower = -1, ifbd = 0, i_fbd = 0; - double fbound = 1.0, f_bounds = 1.0, f_delta_bounds = 1.0; - double ff, y_new, ff_alt; - - for (i = 0; i < neq_; i++) { - y_new = y[i] + step0[i]; - /* - * Force the step to only take 80% a step towards the lower bounds - */ - if (step0[i] < 0.0) { - if (y_new < m_y_low_bounds[i]) { - double legalDelta = 0.8*(m_y_low_bounds[i] - y[i]); - ff = legalDelta / step0[i]; - if (ff < f_bounds) { - f_bounds = ff; - i_lower = i; - } - } - } - /* - * Force the step to only take 80% a step towards the high bounds - */ - if (step0[i] > 0.0) { - if (y_new > m_y_high_bounds[i]) { - double legalDelta = 0.8*(m_y_high_bounds[i] - y[i]); - ff = legalDelta / step0[i]; - if (ff < f_bounds) { - f_bounds = ff; - i_lower = i; - } - } - } - /** - * Now do a delta bounds - * Increase variables by a factor of 2 only - * decrease variables by a factor of 5 only - */ - ff = 1.0; - if ((fabs(y_new) > 2.0 * fabs(y[i])) && - (fabs(y_new-y[i]) > m_ewt[i])) { - ff = fabs(y[i]/(y_new - y[i])); - ff_alt = fabs(m_ewt[i] / (y_new - y[i])); - ff = MAX(ff, ff_alt); - ifbd = 1; - } - if ((fabs(5.0 * y_new) < fabs(y[i])) && - (fabs(y_new - y[i]) > m_ewt[i])) { - ff = y[i]/(y_new-y[i]) * (1.0 - 5.0)/5.0; - ff_alt = fabs(m_ewt[i] / (y_new - y[i])); - ff = MAX(ff, ff_alt); - ifbd = 0; - } - if (ff < f_delta_bounds) { - f_delta_bounds = ff; - i_fbd = ifbd; - } - f_delta_bounds = MIN(f_delta_bounds, ff); - } - fbound = MIN(f_bounds, f_delta_bounds); - /* - * Report on any corrections - */ - if (loglevel > 1) { - if (fbound != 1.0) { - if (f_bounds < f_delta_bounds) { - printf("\t\tboundStep: Variable %d causing bounds " - "damping of %g\n", - i_lower, f_bounds); - } else { - if (ifbd) { - printf("\t\tboundStep: Decrease of Variable %d causing " - "delta damping of %g\n", - i_fbd, f_delta_bounds); - } else { - printf("\t\tboundStep: Increase of variable %d causing" - "delta damping of %g\n", - i_fbd, f_delta_bounds); - } - } - } - } - //return fbound; - return 1.0; - } - - /************************************************************************** - * - * dampStep(): - * - * On entry, step0 must contain an undamped Newton step to the - * current solution y0. This method attempts to find a damping coefficient - * such that the next undamped step would have a norm smaller than - * that of step0. If successful, the new solution after taking the - * damped step is returned in y1, and the undamped step at y1 is - * returned in step1. - */ - int NonlinearSolver::dampStep(const double time_curr, const double* y0, - const double *ydot0, const double* step0, - double* const y1, double* const ydot1, double* step1, - double& s1, SquareMatrix& jac, - int& loglevel, bool writetitle, - int& num_backtracks) { - - - // Compute the weighted norm of the undamped step size step0 - double s0 = solnErrorNorm(step0); - - // Compute the multiplier to keep all components in bounds - // A value of one indicates that there is no limitation - // on the current step size in the nonlinear method due to - // bounds constraints (either negative values of delta - // bounds constraints. - double fbound = boundStep(y0, step0, loglevel); - - // if fbound is very small, then y0 is already close to the - // boundary and step0 points out of the allowed domain. In - // this case, the Newton algorithm fails, so return an error - // condition. - if (fbound < 1.e-10) { - if (loglevel > 1) printf("\t\t\tdampStep: At limits.\n"); - return -3; - } - - //-------------------------------------------- - // Attempt damped step - //-------------------------------------------- - - // damping coefficient starts at 1.0 - double damp = 1.0; - int j, m; - double ff; - num_backtracks = 0; - for (m = 0; m < NDAMP; m++) { - - ff = fbound*damp; - - // step the solution by the damped step size - /* - * Whenever we update the solution, we must also always - * update the time derivative. - */ - for (j = 0; j < neq_; j++) { - y1[j] = y0[j] + ff * step0[j]; - } - calc_ydot(m_order, y1, ydot1); - - doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, y1, ydot1, step1, loglevel); - - // compute the next undamped step, step1[], that would result - // if y1[] were accepted. - - doNewtonSolve(time_curr, y1, ydot1, step1, jac, loglevel); - - // compute the weighted norm of step1 - s1 = solnErrorNorm(step1); - - // write log information - if (loglevel > 3) { - print_solnDelta_norm_contrib((const double *) step0, - "DeltaSolnTrial", - (const double *) step1, - "DeltaSolnTrialTest", - "dampNewt: Important Entries for " - "Weighted Soln Updates:", - y0, y1, ff, 5); - } - if (loglevel > 1) { - printf("\t\t\tdampNewt: s0 = %g, s1 = %g, fbound = %g," - "damp = %g\n", s0, s1, fbound, damp); - } - - - // if the norm of s1 is less than the norm of s0, then - // accept this damping coefficient. Also accept it if this - // step would result in a converged solution. Otherwise, - // decrease the damping coefficient and try again. - - if (s1 < 1.0E-5 || s1 < s0) { - if (loglevel > 2) { - if (s1 > s0) { - if (s1 > 1.0) { - printf("\t\t\tdampStep: current trial step and damping" - " coefficient accepted because test step < 1\n"); - printf("\t\t\t s1 = %g, s0 = %g\n", s1, s0); - } - } - } - break; - } else { - if (loglevel > 1) { - printf("\t\t\tdampStep: current step rejected: (s1 = %g > " - "s0 = %g)", s1, s0); - if (m < (NDAMP-1)) { - printf(" Decreasing damping factor and retrying"); - } else { - printf(" Giving up!!!"); - } - printf("\n"); - } - } - num_backtracks++; - damp /= DampFactor; - } - - // If a damping coefficient was found, return 1 if the - // solution after stepping by the damped step would represent - // a converged solution, and return 0 otherwise. If no damping - // coefficient could be found, return -2. - if (m < NDAMP) { - if (s1 > 1.0) return 0; - else return 1; - } else { - if (s1 < 0.5 && (s0 < 0.5)) return 1; - if (s1 < 1.0) return 0; - return -2; - } - } - - /** - * - * solve_nonlinear_problem(): - * - * Find the solution to F(X) = 0 by damped Newton iteration. On - * entry, x0 contains an initial estimate of the solution. On - * successful return, x1 contains the converged solution. - * - * SolnType = TRANSIENT -> we will assume we are relaxing a transient - * equation system for now. Will make it more general later, - * if an application comes up. - * - */ - int NonlinearSolver::solve_nonlinear_problem(int SolnType, double* y_comm, - double* ydot_comm, double CJ, - double time_curr, - SquareMatrix& jac, - int &num_newt_its, - int &num_linear_solves, - int &num_backtracks, - int loglevelInput) - { - clockWC wc; - - int m = 0; - bool forceNewJac = false; - double s1=1.e30; - - std::vector y_curr(neq_, 0.0); - std::vector ydot_curr(neq_, 0.0); - std::vector stp(neq_, 0.0); - std::vector stp1(neq_, 0.0); - - std::vector y_new(neq_, 0.0); - std::vector ydot_new(neq_, 0.0); - - mdp::mdp_copy_dbl_1(DATA_PTR(y_curr), y_comm, neq_); - // copyn((size_t)neq_, y_comm, y_curr); - mdp::mdp_copy_dbl_1(DATA_PTR(ydot_curr), ydot_comm, neq_); - - - - bool frst = true; - num_newt_its = 0; - num_linear_solves = - m_numTotalLinearSolves; - num_backtracks = 0; - int i_backtracks; - int loglevel = loglevelInput; - - while (1 > 0) { - - /* - * Increment Newton Solve counter - */ - m_numTotalNewtIts++; - num_newt_its++; - - - if (loglevel > 1) { - printf("\t\tSolve_Nonlinear_Problem: iteration %d:\n", - num_newt_its); - } - - // Check whether the Jacobian should be re-evaluated. - - forceNewJac = true; - - if (forceNewJac) { - if (loglevel > 1) { - printf("\t\t\tGetting a new Jacobian and solving system\n"); - } - beuler_jac(jac, DATA_PTR(m_resid), time_curr, CJ, DATA_PTR(y_curr), DATA_PTR(ydot_curr), - num_newt_its); - } else { - if (loglevel > 1) { - printf("\t\t\tSolving system with old jacobian\n"); - } - } - /* - * Go get new scales - */ - setColumnScales(); - - - doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, - DATA_PTR(y_curr), DATA_PTR(ydot_curr), DATA_PTR(stp), loglevel); - - // compute the undamped Newton step - doNewtonSolve(time_curr, DATA_PTR(y_curr), DATA_PTR(ydot_curr), DATA_PTR(stp), - jac, loglevel); - - // damp the Newton step - m = dampStep(time_curr, DATA_PTR(y_curr), DATA_PTR(ydot_curr), - DATA_PTR(stp), DATA_PTR(y_new), DATA_PTR(ydot_new), - DATA_PTR(stp1), s1, jac, loglevel, frst, i_backtracks); - frst = false; - num_backtracks += i_backtracks; - - /* - * Impose the minimum number of newton iterations critera - */ - if (num_newt_its < m_min_newt_its) { - if (m == 1) m = 0; - } - /* - * Impose max newton iteration - */ - if (num_newt_its > 20) { - m = -1; - if (loglevel > 1) { - printf("\t\t\tDampnewton unsuccessful (max newts exceeded) sfinal = %g\n", s1); - } - } - - if (loglevel > 1) { - if (m == 1) { - printf("\t\t\tDampNewton iteration successful, nonlin " - "converged sfinal = %g\n", s1); - } else if (m == 0) { - printf("\t\t\tDampNewton iteration successful, get new" - "direction, sfinal = %g\n", s1); - } else { - printf("\t\t\tDampnewton unsuccessful sfinal = %g\n", s1); - } - } - - // If we are converged, then let's use the best solution possible - // for an end result. We did a resolve in dampStep(). Let's update - // the solution to reflect that. - // HKM 5/16 -> Took this out, since if the last step was a - // damped step, then adding stp1[j] is undamped, and - // may lead to oscillations. It kind of defeats the - // purpose of dampStep() anyway. - // if (m == 1) { - // for (int j = 0; j < neq_; j++) { - // y_new[j] += stp1[j]; - // HKM setting intermediate y's to zero was a tossup. - // slightly different, equivalent results - // #ifdef DEBUG_HKM - // y_new[j] = MAX(0.0, y_new[j]); - // #endif - // } - // } - - bool m_filterIntermediate = false; - if (m_filterIntermediate) { - if (m == 0) { - (void) filterNewStep(time_n, DATA_PTR(y_new), DATA_PTR(ydot_new)); - } - } - // Exchange new for curr solutions - if (m == 0 || m == 1) { - mdp::mdp_copy_dbl_1(DATA_PTR(y_curr), DATA_PTR(y_new), neq_); - calc_ydot(m_order, DATA_PTR(y_curr), DATA_PTR(ydot_curr)); - } - - // convergence - if (m == 1) goto done; - - // If dampStep fails, first try a new Jacobian if an old - // one was being used. If it was a new Jacobian, then - // return -1 to signify failure. - else if (m < 0) { - goto done; - } - } - - done: - mdp::mdp_copy_dbl_1(y_comm, DATA_PTR(y_curr), neq_); - mdp::mdp_copy_dbl_1(ydot_comm, DATA_PTR(ydot_curr), neq_); - - - num_linear_solves += m_numTotalLinearSolves; - - double time_elapsed = wc.secondsWC(); - if (loglevel > 1) { - if (m == 1) { - printf("\t\tNonlinear problem solved successfully in " - "%d its, time elapsed = %g sec\n", - num_newt_its, time_elapsed); - } - } - return m; - } - - /***************************************************************8 - * - * - */ - void NonlinearSolver:: - print_solnDelta_norm_contrib(const double * const solnDelta0, - const char * const s0, - const double * const solnDelta1, - const char * const s1, - const char * const title, - const double * const y0, - const double * const y1, - double damp, - int num_entries) { - int i, j, jnum; - bool used; - double dmax0, dmax1, error, rel_norm; - printf("\t\t%s currentDamp = %g\n", title, damp); - printf("\t\t I ysoln %10s ysolnTrial " - "%10s weight relSoln0 relSoln1\n", s0, s1); - int *imax = mdp::mdp_alloc_int_1(num_entries, -1); - printf("\t\t "); print_line("-", 90); - for (jnum = 0; jnum < num_entries; jnum++) { - dmax1 = -1.0; - for (i = 0; i < neq_; i++) { - used = false; - for (j = 0; j < jnum; j++) { - if (imax[j] == i) used = true; - } - if (!used) { - error = solnDelta0[i] / m_ewt[i]; - rel_norm = sqrt(error * error); - error = solnDelta1[i] / m_ewt[i]; - rel_norm += sqrt(error * error); - if (rel_norm > dmax1) { - imax[jnum] = i; - dmax1 = rel_norm; - } - } - } - if (imax[jnum] >= 0) { - i = imax[jnum]; - error = solnDelta0[i] / m_ewt[i]; - dmax0 = sqrt(error * error); - error = solnDelta1[i] / m_ewt[i]; - dmax1 = sqrt(error * error); - printf("\t\t %4d %12.4e %12.4e %12.4e %12.4e " - "%12.4e %12.4e %12.4e\n", - i, y0[i], solnDelta0[i], y1[i], - solnDelta1[i], m_ewt[i], dmax0, dmax1); - } - } - printf("\t\t "); print_line("-", 90); - mdp::mdp_safe_free((void **) &imax); - } - - - -} - diff --git a/Cantera/src/numerics/NonlinearSolver.h b/Cantera/src/numerics/NonlinearSolver.h deleted file mode 100644 index ef2d13c06..000000000 --- a/Cantera/src/numerics/NonlinearSolver.h +++ /dev/null @@ -1,329 +0,0 @@ -/** - * @file NonlinearSolve.h - * Class that calculates the solution to a nonlinear, dense, set - * of equations (see \ref numerics - * and class \link Cantera::NonlinearSolver NonlinearSolver\endlink). - */ -/* - * Copywrite 2004 Sandia Corporation. Under the terms of Contract - * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government - * retains certain rights in this software. - * See file License.txt for licensing information. - */ - -#ifndef CT_NONLINEARSOLVER_H -#define CT_NONLINEARSOLVER_H - -#include "ResidJacEval.h" - -namespace Cantera { - -#define NSOLN_TYPE_PSEUDO_TIME_DEPENDENT 2 -#define NSOLN_TYPE_TIME_DEPENDENT 1 -#define NSOLN_TYPE_STEADY_STATE 0 - - //! Class that calculates the solution to a nonlinear system - /*! - * - * @ingroup numerics - */ - class NonlinearSolver { - - //! Default constructor - /*! - * @param func Residual and jacobian evaluator function object - */ - NonlinearSolver(ResidJacEval *func); - - //!Copy Constructor for the %ThermoPhase object. - /*! - * @param right Item to be copied - */ - NonlinearSolver(const NonlinearSolver &right); - - //! Destructor - ~NonlinearSolver(); - - //! Assignment operator - /*! - * This is NOT a virtual function. - * - * @param right Reference to %NonlinearSolver object to be - * copied into the - * current one. - */ - NonlinearSolver& operator=(const NonlinearSolver &right); - - - //! Create solution weights for convergence criteria - /*! - * We create soln weights from the following formula - * - * wt[i] = rtol * abs(y[i]) + atol[i] - * - * The program always assumes that atol is specific - * to the solution component - * - * param y vector of the current solution values - */ - void createSolnWeights(const double * const y); - - - //! L2 norm of the delta of the solution vector - /*! - * calculate the norm of the solution vector. This will - * involve the column scaling of the matrix - * - * The second argument has a default of false. However, - * if true, then a table of the largest values is printed - * out to standard output. - */ - double solnErrorNorm(const double * const delta_y, - bool printLargest = false); - - //! L2 norm of the residual of the equation system - /*! - * Calculate the norm of the residual vector. This may - * involve using the row sum scaling from the matrix problem. - * - * The second argument has a default of false. However, - * if true, then a table of the largest values is printed - * out to standard output. - */ - double residErrorNorm(const double * const resid, - bool printLargest = false); - - //! Compute the current Residual - /*! - * Compute the time dependent residual of - * the set of equations. - */ - void doTDResidualCalc(const double time_curr, const int typeCalc, - const double * const y_curr, - const double * const ydot_curr, double* const residual, - int loglevel); - - //! Compute the current Residual - /*! - * Compute the steady state residual of - * the set of equations. - */ - void doSteadyResidualCalc(const double time_curr, const int typeCalc, - const double * const y_curr, - double* const residual, int loglevel); - - void doResidualCalc(const double time_curr, const int typeCalc, - const double * const y_curr, - const double * const ydot_curr, double* const residual, - int loglevel); - - //! Compute the undamped Newton step - /*! - * - * Compute the undamped Newton step. The residual function is - * evaluated at the current time, t_n, at the current values of the - * solution vector, m_y_n, and the solution time derivative, m_ydot_n. - * The Jacobian is not recomputed. - * - * A factored jacobian is reused, if available. If a factored jacobian - * is not available, then the jacobian is factored. Before factoring, - * the jacobian is row and column-scaled. Column scaling is not - * recomputed. The row scales are recomputed here, after column - * scaling has been implemented. - * - */ - void doNewtonSolve(const double time_curr, const double * const y_curr, - const double * const ydot_curr, double* const delta_y, - SquareMatrix& jac, int loglevel); - - //! - /*! - * - * Return the factor by which the undamped Newton step 'step0' - * must be multiplied in order to keep all solution components in - * all domains between their specified lower and upper bounds. - * Other bounds may be applied here as well. - * - * Currently the bounds are hard coded into this routine: - * - * Minimum value for all variables: - 0.01 * m_ewt[i] - * Maximum value = none. - * - * Thus, this means that all solution components are expected - * to be numerical greater than zero in the limit of time step - * truncation errors going to zero. - * - * Delta bounds: The idea behind these is that the Jacobian - * couldn't possibly be representative if the - * variable is changed by a lot. (true for - * nonlinear systems, false for linear systems) - * Maximum increase in variable in any one newton iteration: - * factor of 2 - * Maximum decrease in variable in any one newton iteration: - * factor of 5 - */ - double boundStep(const double* const y, - const double* const step0, const int loglevel); - - - //! set bounds constraints for all variables in the problem - /*! - * - * @param y_low_bounds Vector of lower bounds - * @param y_high_bounds Vector of high bounds - */ - void setBoundsConstraints(const double * const y_low_bounds, - const double * const y_high_bounds); - - /** - * Internal function to calculate the predicted solution - * at a time step. - */ - void calc_y_pred(int); - - /** - * Internal function to calculate the time derivative at the - * new step - */ - void calc_ydot(int order, double * const y_curr, double * const ydot_curr); - - void beuler_jac(SquareMatrix &, double * const, - double, double, double * const, double * const, int); - - - double filterNewStep(double, double *, double *); - - //! Find a damping coefficient through a look-ahead mechanism - /*! - * On entry, step0 must contain an undamped Newton step for the - * solution x0. This method attempts to find a damping coefficient - * such that all components stay in bounds, and the next - * undamped step would have a norm smaller than - * that of step0. If successful, the new solution after taking the - * damped step is returned in y1, and the undamped step at y1 is - * returned in step1. - * - * @param time_curr Current physical time - * @param y0 Base value of the solution before any steps - * are taken - * @param ydot0 Base value of the time derivative of teh - * solution - * @param step0 Initial step suggested. - * @param y1 - */ - int dampStep(const double time_curr, const double* y0, - const double *ydot0, const double* step0, - double* const y1, double* const ydot1, double* step1, - double& s1, SquareMatrix& jac, - int& loglevel, bool writetitle, - int& num_backtracks); - - - - // Compute the weighted norm of the undamped step size step0 - - //! Find the solution to F(X) = 0 by damped Newton iteration. - /*! - * On - * entry, x0 contains an initial estimate of the solution. On - * successful return, x1 contains the converged solution. - * - * SolnType = TRANSIENT -> we will assume we are relaxing a transient - * equation system for now. Will make it more general later, - * if an application comes up. - * - */ - int solve_nonlinear_problem(int SolnType, double* y_comm, - double* ydot_comm, double CJ, - double time_curr, - SquareMatrix& jac, - int &num_newt_its, - int &num_linear_solves, - int &num_backtracks, - int loglevelInput); - - - void setColumnScales(); - - void - print_solnDelta_norm_contrib(const double * const solnDelta0, - const char * const s0, - const double * const solnDelta1, - const char * const s1, - const char * const title, - const double * const y0, - const double * const y1, - double damp, - int num_entries); - - - - //! Pointer to the residual and jacobian evaluator for the - //! function - /*! - * See ResidJacEval.h for an evaluator. - */ - ResidJacEval *m_func; - - //! Local copy of the number of equations - int neq_; - - //! Soln error weights - std::vector m_ewt; - - std::vector m_y_n; - std::vector m_y_nm1; - std::vector m_colScales; - - //! Weights for normalizing the values of the residuals - - std::vector m_rowScales; - - std::vector m_resid; - - //! Bounds vector for each species - std::vector m_y_high_bounds; - - //! Lower bounds vector for each species - std::vector m_y_low_bounds; - - double delta_t_n; - - //! Counter for the total number of function evaluations - int m_nfe; - - //! The type of column scaled used in the solution of the problem - bool m_colScaling; - - //! int indicating whether row scaling is turned on (1) or not (0) - int m_rowScaling; - - int m_numTotalLinearSolves; - - int m_numTotalNewtIts; - - int m_min_newt_its; - - int filterNewstep; - - //! Current system time - /*! - * Note, we assume even for steady state problems that the residual - * is a function of a system time. - */ - double time_n; - - int m_matrixConditioning; - - int m_order; - - doublereal rtol_; - - doublereal atolBase_; - - std::vector atolk_; - }; - -} - -#endif diff --git a/Cantera/src/numerics/ResidEval.h b/Cantera/src/numerics/ResidEval.h deleted file mode 100755 index 9f8315649..000000000 --- a/Cantera/src/numerics/ResidEval.h +++ /dev/null @@ -1,102 +0,0 @@ -/** - * @file ResidEval.h - * - */ -// Copyright 2006 California Institute of Technology - -#ifndef CT_RESIDEVAL_H -#define CT_RESIDEVAL_H - -#include "ct_defs.h" -#include "ctexceptions.h" - -namespace Cantera { - - const int c_NONE = 0; - const int c_GE_ZERO = 1; - const int c_GT_ZERO = 2; - const int c_LE_ZERO = -1; - const int c_LT_ZERO = -2; - - /** - * Virtual base class for DAE residual function evaluators. - * Classes derived from ResidEval evaluate the residual function - * \f[ - \vec{F}(t,\vec{y}, \vec{y^\prime}) - * \f] - * The DAE solver attempts to find a solution y(t) such that F = 0. - * @ingroup DAE_Group - */ - class ResidEval { - - public: - - ResidEval() {} - virtual ~ResidEval() {} - - /** - * Constrain solution component k. Possible values for - * 'flag' are: - * - c_NONE no constraint - * - c_GE_ZERO >= 0 - * - c_GT_ZERO > 0 - * - c_LE_ZERO <= 0 - * - c_LT_ZERO < 0 - */ - virtual void constrain(const int k, const int flag) { m_constrain[k] = flag; } - int constraint(const int k) const { - std::map::const_iterator i = m_constrain.find(k); - if (i != m_constrain.end()) return i->second; - return c_NONE; - } - - /** - * Specify that solution component k is purely algebraic - - * that is, the derivative of this component does not appear - * in the residual function. - */ - virtual void setAlgebraic(const int k) { m_alg[k] = 1; } - virtual bool isAlgebraic(const int k) {return (m_alg[k] == 1); } - - - /** - * Evaluate the residual function. Called by the - * integrator. - * @param t time. (input) - * @param y solution vector. (input) - * @param ydot rate of change of solution vector. (input) - * @param r residual vector (output) - */ - virtual int eval(const doublereal t, const doublereal * const y, - const doublereal * const ydot, - doublereal * const r) { - throw CanteraError("ResidEval::eval()", "base class called"); - } - - /** - * Fill the solution and derivative vectors with the initial - * conditions at initial time t0. If these do not satisfy the - * residual equation, call one of the "corrrectInitial_xxx" - * methods before calling solve. - */ - virtual void getInitialConditions(const doublereal t0, doublereal * const y, - doublereal * const ydot) { - throw CanteraError("ResidEval::GetInitialConditions()", "base class called"); - } - - //! Return the number of equations in the equation system - virtual int nEquations() const = 0; - - - protected: - - std::map m_alg; - std::map m_constrain; - - private: - - }; - -} - -#endif diff --git a/Cantera/src/numerics/ResidJacEval.cpp b/Cantera/src/numerics/ResidJacEval.cpp deleted file mode 100644 index 97f8bb5ed..000000000 --- a/Cantera/src/numerics/ResidJacEval.cpp +++ /dev/null @@ -1,282 +0,0 @@ -/** - * @file ResidJacEval.cpp - * - */ -/* - * Copywrite 2004 Sandia Corporation. Under the terms of Contract - * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government - * retains certain rights in this software. - * See file License.txt for licensing information. - */ - -#include "ct_defs.h" -#include "ctlapack.h" -#include "ResidJacEval.h" - -#include -#include - -using namespace std; - -namespace Cantera { - - /************************************************************************* - * - * ResidJacEval(): - * - * Default constructor for the ResidJacEval class. - * - * atol has a default of 1.0E-13. - */ - ResidJacEval::ResidJacEval(doublereal atol) : - ResidEval(), - m_atol(atol) - { - } - - // Copy Constructor for the %ResidJacEval object - /* - */ - ResidJacEval::ResidJacEval(const ResidJacEval &right) : - ResidEval() - { - *this = operator=(right); - } - - /* - * - */ - ResidJacEval::~ResidJacEval() - { - } - - ResidJacEval& ResidJacEval::operator=(const ResidJacEval &right) { - if (this == &right) { - return *this; - } - - ResidEval::operator=(right); - - m_atol = right.m_atol; - neq_ = right.neq_; - - return *this; - } - - // Duplication routine for objects which inherit from - // %ResidJacEval - /* - * This virtual routine can be used to duplicate %ResidJacEval objects - * inherited from %ResidJacEval even if the application only has - * a pointer to %ResidJacEval to work with. - * - * These routines are basically wrappers around the derived copy - * constructor. - */ - ResidJacEval *ResidJacEval::duplMyselfAsResidJacEval() const { - ResidJacEval *ff = new ResidJacEval(*this); - return ff; - } - - int ResidJacEval::nEquations() const { - return neq_; - } - - /* - * - * setAtol(): - * - * Set the absolute tolerance value - */ - void ResidJacEval::setAtol(doublereal atol) - { - m_atol = atol; - if (m_atol <= 0.0) { - throw CanteraError("ResidJacEval::setAtol", - "atol must be greater than zero"); - } - } - - /************************************************************************** - * - * - * - * Fill the solution vector with the initial conditions - * at initial time t0. - */ - void ResidJacEval:: - getInitialConditionsDot(const doublereal t0, const size_t leny, - doublereal * const y, doublereal * const ydot) { - for (int i = 0; i < neq_; i++) { - y[i] = 0.0; - } - if (ydot) { - for (int i = 0; i < neq_; i++) { - ydot[i] = 0.0; - } - } - } - - /************************************************************************** - * - * - * - * Fill the solution vector with the initial conditions - * at initial time t0. - * - */ - void ResidJacEval:: - getInitialConditions(doublereal t0, - doublereal * const y, doublereal * const ydot) { - size_t leny = neq_; - getInitialConditionsDot(t0, leny, y, 0); - } - - /************************************************************************** - * - * user_out(): - * - * This function may be used to create output at various points in the - * execution of an application. - * - */ - void ResidJacEval:: - user_out2(const int ifunc, const doublereal t, const doublereal deltaT, - const doublereal *y, const doublereal *ydot) { - - } - - void ResidJacEval:: - user_out(const int ifunc, const doublereal t, - const doublereal *y, const doublereal *ydot) { - user_out2(ifunc, t, 0.0, y, ydot); - } - - /************************************************************************** - * - * - */ - void ResidJacEval:: - evalTimeTrackingEqns(const doublereal t, const doublereal deltaT, - const doublereal *y, - const doublereal *ydot) { - - } - - /******************************************************************** - * - * - * - * Return a vector of delta y's for calculation of the - * numerical Jacobian - */ - void ResidJacEval:: - calcDeltaSolnVariables(const doublereal t, - const doublereal * const ySoln, - const doublereal * const ySolnDot, - doublereal * const deltaYSoln, - const doublereal *const solnWeights) - { - if (!solnWeights) { - for (int i = 0; i < neq_; i++) { - deltaYSoln[i] = m_atol + fabs(1.0E-6 * ySoln[i]); - } - } else { - for (int i = 0; i < neq_; i++) { - deltaYSoln[i] = m_atol + - fmaxx(1.0E-2 * solnWeights[i], 1.0E-6 * fabs(ySoln[i])); - } - } - } - - /****************************************************************** - * - * calcSolnScales(): - * - * Returns a vector of ysolnScales[] that can be used to column scale - * Jacobians. - */ - void ResidJacEval:: - calcSolnScales(const doublereal t, - const doublereal * const ysoln, - const doublereal * const ysolnOld, - doublereal * const ysolnScales) - { - for (int i = 0; i < neq_; i++) { - ysolnScales[i] = 1.0; - } - } - - void ResidJacEval::filterSolnPrediction(doublereal t, - doublereal * const y) { - - } - - /************************************************************************** - * - * evalStoppingCriteria() - * - * If there is a stopping critera other than time set it here. - * - */ - bool ResidJacEval:: - evalStoppingCritera(doublereal &time_current, - doublereal &delta_t_n, - doublereal *y_n, - doublereal *ydot_n) - { - return false; - } - - /************************************************************************** - * - * matrixConditioning() - * - * Multiply the matrix by the inverse of a matrix which lead to a - * better conditioned system. The default, specified here, is to - * do nothing. - */ - void ResidJacEval:: - matrixConditioning(doublereal * const matrix, const int nrows, - doublereal * const rhs) - { - } - - /************************************************************************** - * - */ - void ResidJacEval:: - evalResidNJ(doublereal t, const doublereal deltaT, - const doublereal * y, - const doublereal * ydot, - doublereal * resid, - bool NJevaluation, - int id_x, - doublereal delta_x) - { - printf("Not implemented\n"); - std::exit(-1); - } - - /************************************************************************** - * - * evalJacobian() - * - * Calculate the jacobian and the residual at the current - * time and values. - * Backwards Euler is assumed. - */ - void ResidJacEval:: - evalJacobian(const doublereal t, const doublereal deltaT, - const doublereal * const y, - const doublereal * const ydot, - SquareMatrix &J, - doublereal * const resid) - { - printf("Not implemented\n"); - std::exit(-1); - } - - -} - diff --git a/Cantera/src/numerics/ResidJacEval.h b/Cantera/src/numerics/ResidJacEval.h deleted file mode 100644 index fd9a89fc9..000000000 --- a/Cantera/src/numerics/ResidJacEval.h +++ /dev/null @@ -1,175 +0,0 @@ -/** - * @file ResidJacEval.h - * - * Dense, Square (not sparse) matrices. - */ -/* - * Copywrite 2004 Sandia Corporation. Under the terms of Contract - * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government - * retains certain rights in this software. - * See file License.txt for licensing information. - */ - - -#ifndef CT_RESIDJACEVAL_H -#define CT_RESIDJACEVAL_H - -#include "ResidEval.h" -#include "SquareMatrix.h" - -namespace Cantera { - - /** - * A class for full (non-sparse) matrices with Fortran-compatible - * data storage. Adds matrix operations to class Array2D. - */ - class ResidJacEval : public ResidEval { - - public: - - /** - * Default constructor - */ - ResidJacEval(doublereal atol = 1.0e-13); - - //!Copy Constructor for the %ResidJacEval object - /*! - * @param right Item to be copied - */ - ResidJacEval(const ResidJacEval &right); - - /// Destructor. Does nothing. - virtual ~ResidJacEval(); - - //! Assignment operator - /*! - * This is NOT a virtual function. - * - * @param right Reference to %ResidJacEval object to be copied into the - * current one. - */ - ResidJacEval& operator=(const ResidJacEval &right); - - //! Duplication routine for objects which inherit from - //! residJacEval - /*! - * This virtual routine can be used to duplicate %ResidJacEval objects - * inherited from %ResidJacEval even if the application only has - * a pointer to %ResidJacEval to work with. - * - * These routines are basically wrappers around the derived copy - * constructor. - */ - virtual ResidJacEval *duplMyselfAsResidJacEval() const; - - //! Return the number of equations in the equation system - virtual int nEquations() const; - - /** - * Evaluate the residual function. - * @param t time (input, do not modify) - * @param y solution vector (input, do not modify) - * @param ydot rate of change of solution vector. (input, do - * not modify) - */ - virtual void evalResidNJ(doublereal t, const doublereal deltaT, - const doublereal * const y, - const doublereal * const ydot, - doublereal * const resid, - bool NJevaluation = false, - int id_x = 0, - doublereal delta_x = 0.0); - - /** - * Fill the solution vector with the initial conditions - * at initial time t0. - */ - virtual void getInitialConditionsDot(const doublereal t0, size_t leny, - doublereal * const y, - doublereal * const ydot); - - virtual void getInitialConditions(const doublereal t0, - doublereal * const y, - doublereal * const ydot); - - virtual void filterSolnPrediction(doublereal t, - doublereal * const y); - - void setAtol(doublereal atol); - - virtual void evalTimeTrackingEqns(const doublereal t, const doublereal deltaT, - const doublereal * const y, - const doublereal * const ydot); - - virtual bool evalStoppingCritera(doublereal &time_current, - doublereal &delta_t_n, - doublereal *y_n, - doublereal *ydot_n); - /** - * Return a vector of delta y's for calculation of the - * numerical Jacobian - */ - virtual void - calcDeltaSolnVariables(const doublereal t, - const doublereal * const ysoln, - const doublereal * const ysolnDot, - doublereal * const deltaYsoln, - const doublereal * const solnWeights=0); - - /** - * Returns a vector of ysolnScales[] that can be used to column - * scale Jacobians. - */ - virtual void calcSolnScales(const doublereal t, - const doublereal * const ysoln, - const doublereal * const ysolnOld, - doublereal * const ysolnScales); - - /** - * This function may be used to create output at various points in the - * execution of an application. - * - */ - virtual void user_out2(const int ifunc, const doublereal t, - const doublereal deltaT, - const doublereal * const y, - const doublereal * const ydot); - - virtual void user_out(const int ifunc, const doublereal t, - const doublereal *y, - const doublereal *ydot); - - - virtual void matrixConditioning(doublereal * const matrix, const int nrows, - doublereal * const rhs); - - /********************************************************************* - * - * evalJacobian() - * - * Calculate the jacobian and the residual at the current - * time and values. - * Backwards Euler is assumed. - */ - virtual void evalJacobian(const doublereal t, const doublereal deltaT, - - const double* const y, - const double* const ydot, - SquareMatrix &J, - doublereal * const resid); - - - - protected: - - doublereal m_atol; - - int neq_; - - }; -} - -#endif - - - diff --git a/Cantera/src/numerics/SquareMatrix.cpp b/Cantera/src/numerics/SquareMatrix.cpp deleted file mode 100644 index ccf7749ae..000000000 --- a/Cantera/src/numerics/SquareMatrix.cpp +++ /dev/null @@ -1,114 +0,0 @@ -/** - * @file DenseMatrix.cpp - * - */ -/* - * Copywrite 2004 Sandia Corporation. Under the terms of Contract - * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government - * retains certain rights in this software. - * See file License.txt for licensing information. - */ - -#include "ct_defs.h" -#include "stringUtils.h" -#include "ctlapack.h" -#include "SquareMatrix.h" - -#include -#include - -#include - -using namespace std; - -namespace Cantera { - /** - * - * copy constructor - */ - SquareMatrix::SquareMatrix(const SquareMatrix& y) : - DenseMatrix(y), - m_factored(y.m_factored) - { - } - - /** - * Assignment operator - */ - SquareMatrix& SquareMatrix::operator=(const SquareMatrix& y) { - if (&y == this) return *this; - DenseMatrix::operator=(y); - m_factored = y.m_factored; - return *this; - } - - /** - * Solve Ax = b. Vector b is overwritten on exit with x. - */ - int SquareMatrix::solve(double* b) - { - int info=0; - /* - * Check to see whether the matrix has been factored. - */ - if (!m_factored) { - factor(); - } - /* - * Solve the factored system - */ - ct_dgetrs(ctlapack::NoTranspose, static_cast(nRows()), - 1, &(*(begin())), static_cast(nRows()), - DATA_PTR(ipiv()), b, static_cast(nColumns()), info); - if (info != 0) - throw CanteraError("SquareMatrix::solve", - "DGETRS returned INFO = "+int2str(info)); - return 0; - } - - /** - * Set all entries to zero - */ - void SquareMatrix::zero() { - int n = static_cast(nRows()); - if (n > 0) { - int nn = n * n; - double *sm = &m_data[0]; - /* - * Using memset is the fastest way to zero a contiguous - * section of memory. - */ - (void) memset((void *) sm, 0, nn * sizeof(double)); - } - } - - /** - * Factor A. A is overwritten with the LU decomposition of A. - */ - int SquareMatrix::factor() { - integer n = static_cast(nRows()); - int info=0; - m_factored = true; - ct_dgetrf(n, n, &(*(begin())), static_cast(nRows()), - DATA_PTR(ipiv()), info); - if (info != 0) { - cout << "Singular matrix, info = " << info << endl; - throw CanteraError("invert", - "DGETRF returned INFO="+int2str(info)); - } - return 0; - } - /* - * clear the factored flag - */ - void SquareMatrix::clearFactorFlag() { - m_factored = false; - } - /** - * set the factored flag - */ - void SquareMatrix::setFactorFlag() { - m_factored = true; - } -} - diff --git a/Cantera/src/numerics/SquareMatrix.h b/Cantera/src/numerics/SquareMatrix.h deleted file mode 100644 index d549d8781..000000000 --- a/Cantera/src/numerics/SquareMatrix.h +++ /dev/null @@ -1,91 +0,0 @@ -/** - * @file SquareMatrix.h - * Dense, Square (not sparse) matrices. - */ -/* - * Copywrite 2004 Sandia Corporation. Under the terms of Contract - * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government - * retains certain rights in this software. - * See file License.txt for licensing information. - */ - -#ifndef CT_SQUAREMATRIX_H -#define CT_SQUAREMATRIX_H - -#include "DenseMatrix.h" - -namespace Cantera { - - /** - * A class for full (non-sparse) matrices with Fortran-compatible - * data storage. Adds matrix inversion operations to this class from DenseMatrix. - */ - class SquareMatrix: public DenseMatrix { - - public: - - SquareMatrix(): - DenseMatrix(), - m_factored(false) - { - } - - /** - * Constructor. Create an \c n by \c n matrix, and initialize - * all elements to \c v. - */ - SquareMatrix(int n, doublereal v = 0.0) : - DenseMatrix(n, n, v), - m_factored(false) - { - } - - /** - * Copy Constructor - */ - SquareMatrix(const SquareMatrix&); - - /** - * Assignment operator - */ - SquareMatrix& operator=(const SquareMatrix&); - - - /// Destructor. Does nothing. - virtual ~SquareMatrix(){} - - /** - * Solves the Ax = b system returning x in the b spot. - */ - int solve(double *b); - - /** - * Zero the matrix - */ - void zero(); - - /** - * Factors the A matrix, overwriting A. We flip m_factored - * boolean to indicate that the matrix is now A-1. - */ - int factor(); - /** - * clear the factored flag - */ - void clearFactorFlag(); - /** - * set the factored flag - */ - void setFactorFlag(); - - /* - * the factor flag - */ - bool m_factored; - }; -} - -#endif - - -