research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Removed some unused numerical integration utilities

Browse source
This commit is contained in:
Ray Speth 2012-01-18 23:44:50 +00:00
parent d17ba8fbcb
commit 6589eb7f7c
9 changed files with 0 additions and 4574 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1784
Cantera/src/base/mdp_allo.cpp
View file

File diff suppressed because it is too large Load diff

734
Cantera/src/base/mdp_allo.h
View file

@ -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 <stdexcept>
/*
* Include the header here in order to pick up size_t definition
*/
#include <cstring>
/**
* 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
/****************************************************************************/

963
Cantera/src/numerics/NonlinearSolver.cpp
View file

@ -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 <limits>
#include "SquareMatrix.h"
#include "NonlinearSolver.h"
#include "clockWC.h"
#include "vec_functions.h"
#include <ctime>
#include "mdp_allo.h"
#include <cfloat>
extern void print_line(const char *, int);
#include <vector>
#include <cstdio>
#include <cmath>
#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<double>::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<double> 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<doublereal> y_curr(neq_, 0.0);
std::vector<doublereal> ydot_curr(neq_, 0.0);
std::vector<doublereal> stp(neq_, 0.0);
std::vector<doublereal> stp1(neq_, 0.0);
std::vector<doublereal> y_new(neq_, 0.0);
std::vector<doublereal> 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);
}
}

329
Cantera/src/numerics/NonlinearSolver.h
View file

@ -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<doublereal> m_ewt;
std::vector<doublereal> m_y_n;
std::vector<doublereal> m_y_nm1;
std::vector<doublereal> m_colScales;
//! Weights for normalizing the values of the residuals
std::vector<doublereal> m_rowScales;
std::vector<doublereal> m_resid;
//! Bounds vector for each species
std::vector<doublereal> m_y_high_bounds;
//! Lower bounds vector for each species
std::vector<doublereal> 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<doublereal> atolk_;
};
}
#endif

102
Cantera/src/numerics/ResidEval.h
View file

@ -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<int,int>::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<int, int> m_alg;
std::map<int, int> m_constrain;
private:
};
}
#endif

282
Cantera/src/numerics/ResidJacEval.cpp
View file

@ -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 <iostream>
#include <vector>
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);
}
}

175
Cantera/src/numerics/ResidJacEval.h
View file

@ -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

114
Cantera/src/numerics/SquareMatrix.cpp
View file

@ -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 <iostream>
#include <vector>
#include <cstring>
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<int>(nRows()),
1, &(*(begin())), static_cast<int>(nRows()),
DATA_PTR(ipiv()), b, static_cast<int>(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<int>(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<int>(nRows());
int info=0;
m_factored = true;
ct_dgetrf(n, n, &(*(begin())), static_cast<int>(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;
}
}

91
Cantera/src/numerics/SquareMatrix.h
View file

@ -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