/*
Jacobian related routines for the NewtonLib package.
* Written by L. Weimann
* Purpose Jacobian and linear algebra related routines
* Category ???. - Utilities
* Keywords Linear system solution, Numerical differentiation,
Jacobian scaling
* Version 1.1
* Revision May 2006
* Latest Change May 2006
* Library NewtonLib
* Code C, Double Precision
* Environment Standard C environment on PC's,
workstations and hosts.
* Copyright (c) Konrad-Zuse-Zentrum fuer
Informationstechnik Berlin (ZIB)
Takustrasse 7, D-14195 Berlin-Dahlem
phone : + 49/30/84185-0
fax : + 49/30/84185-125
* Contact Bodo Erdmann
ZIB, Division Scientific Computing,
Department Numerical Analysis and Modelling
phone : + 49/30/84185-185
fax : + 49/30/84185-107
e-mail: erdmann@zib.de
---------------------------------------------------------------
* Licence
You may use or modify this code for your own non commercial
purposes for an unlimited time.
In any case you should not deliver this code without a special
permission of ZIB.
In case you intend to use the code commercially, we oblige you
to sign an according licence agreement with ZIB.
* Warranty
This code has been tested up to a certain level. Defects and
weaknesses, which may be included in the code, do not establish
any warranties by ZIB. ZIB does not take over any liabilities
which may follow from acquisition or application of this code.
* Software status
This code is under care of ZIB and belongs to ZIB software class 2.
------------------------------------------------------------
The implementation of the routines within this file depend on the
chosen Jacobian format and the selected linear system solver.
The JACOBIAN type used within this file is defined in the header
file "jacobian.h". By changing the contents of this header file,
you may be able to implement any type of Jacobians, for example,
sparse Jacobians. The current implementation supports Jacobians
stored in full storage mode, and uses the Fortran LAPACK solver
routines DGETRF for the LU-decomposition and DGETRS for the
solution (i.e. backsubstitution).
The routines and functions are the following:
---------------------------------------------
int nleq_jacalloc(int n, JACOBIAN **jac, int *ldjac,
struct NLEQ_OPT *opt)
This function must allocate memory for the Jacobian, using either
directly the malloc function or the functions nleq_fwalloc,
nleq_iwalloc, nleq_pfwalloc (see file "utils.c"), as suitable.
The parameters are:
int n (input) The number of equations and unknowns
JACOBIAN **jac (output) A pointer of type (JACOBIAN *) returned
by the memory allocation function, and
pointing to the reserved memory.
int *ldjac (output) A positive int number specifying the
leading dimension size of the Jacobian
(asuming the Jacobian to be a two dimensional
array of (double)).
struct NLEQ_OPT *opt (input) The input options of the main routine.
The int return code of nleq_jacalloc must be 0, if everything worked fine
with the memory allocation, and otherwise a number from -999 to -996.
---
Note: You may expand the struct NLEQ_OPT, which is defined in the "nleq.h"
header file, by adding fields which may describe characteristics
of your Jacobian, for example, such as lower and upper bandwidth
of a Jacobian which will be stored in band mode.
---
int nleq_linfact(int n, JACOBIAN *jac, struct NLEQ_OPT *opt)
This function must call the linear algebra routine for computation of
the (LU-)decomposition.
The parameters are:
int n (input) The number of equations and unknowns
JACOBIAN *jac (in/out) On input: The Jacobian
On output: The decomposition of the Jacobian.
struct NLEQ_OPT *opt (input) The input options of the main routine.
The int return code must be 0, if the matrix decomposition was successful.
A singular matrix must be indicated by a positive return code, any other
error by a negative return code.
---
Note: To keep all needed information of the matrix decomposition, it
may be necessary to allocate additional memory on the first call
of nleq_linfact, i.e. to store pivot indices or other information.
---
int nleq_linsol(int n, double *b, struct NLEQ_OPT *opt)
This function must call the linear algebra routine for the solution
(i.e. backsubstitution) of the linear system, after the Jacobian has
been factorized through nleq_linfact.
The parameters are:
int n (input) The number of equations and unknowns
double *b (in/out) On input: The right hand side vector of the
linear equations system.
On output: The solution of the linear system.
struct NLEQ_OPT *opt (input) The input options of the main routine.
A successful completion must be indicated by a return code 0, a failure
condition by any nonzero return code.
---
void nleq_linfree()
This routine must free any memory, which has been allocated either by
the routine nleq_jacalloc or the routine nleq_linfact.
---
int nleq_numjac(NLEQ_FFUN *f, int n, double *x, double *fx,
double *scale, JACOBIAN *jac, int *nfcn,
struct NLEQ_OPT *opt)
This function should compute the Jacobian by numerical difference
approximation. The main codes will also work if this function just
does nothing more than returning a nonzero value, but for this case
the Jacobian routine must be always supplied by the user.
The parameters are:
NLEQ_FFUN *f (input) A pointer to the user problem function f.
int n (input) The number of equations and unknowns
and size of any vectors.
double *x (input) The input vector where to compute the Jacobian.
The content of x must be preserved on output!
double *fx (input) The vector f(x).
The content of fx must be preserved on output!
double *scale (input) The scaling vector.
The content of scale must be preserved on
output!
JACOBIAN *jac (output) A pointer to the Jacobian storage.
int *nfcn (in/out) The value of *nfcn must be incremented by one
for each call of the user problem function.
struct NLEQ_OPT *opt (input) The input options of the main routine.
A successful completion must be indicated by a return code 0, a failure
condition by any nonzero return code.
---
void nleq_jacrow_scale(int n, JACOBIAN *jac, double *scale,
struct NLEQ_OPT *opt)
This routine must scale the rows of the Jacobian and change the sign
of all Jacobian elements. The function is used by the main codes
nleq_res and qnres.
The scaling must accomplish the following transformation:
jac[i][j] = - jac[i][j]/scale[i] for i=0,...,n-1 and j=0,...,n-1.
The parameters are:
int n (input) The number of equations and unknowns
and size of any vectors.
JACOBIAN *jac (in/out) A pointer to the Jacobian storage.
double *scale (input) The scaling vector.
The content of scale must be preserved on
output!
struct NLEQ_OPT *opt (input) The input options of the main routine.
---
void nleq_jaccolumn_scale(int n, JACOBIAN *jac, double *scale,
struct NLEQ_OPT *opt)
This routine must scale the columns of the Jacobian and change the sign
of all Jacobian elements. The funtion is used by the main codes
nleq_err and qnerr.
The scaling must accomplish the following transformation:
jac[i][j] = - jac[i][j]*scale[j] for j=0,...,n-1 and i=0,...,n-1.
The parameters are:
int n (input) The number of equations and unknowns
and size of any vectors.
JACOBIAN *jac (in/out) A pointer to the Jacobian storage.
double *scale (input) The scaling vector.
The content of scale must be preserved on
output!
struct NLEQ_OPT *opt (input) The input options of the main routine.
*/
#include
#include
#include
#include "nleq.h"
#ifdef FMAT
#define JMODE 'N'
#else
#define JMODE 'T'
#endif
#define AJDEL 1.0e-8
#define AJMIN 1.0e-4
extern void dgetrf_(int *m, int *n, double *a, int *lda,
int *ipiv, int *info);
extern void dgetrs_(char *trans, int *n, int *nrhs, double *a,
int *lda, int *ipiv, double *b,
int *ldb, int *info);
static int *pivot=NULL;
static JACOBIAN *mat=NULL;
static double *v=NULL;
int nleq_jacalloc(int n, JACOBIAN **jac, int *ldjac,
struct NLEQ_OPT *opt)
{ *ldjac = n;
return nleq_fwalloc(n*n,jac,"jac");
}
int nleq_linfact(int n, JACOBIAN *jac, struct NLEQ_OPT *opt)
{ int fail=0;
if (!pivot)
{ fail=nleq_iwalloc(n,&pivot,"pivot");
if ( fail !=0 ) return fail;
};
mat=jac;
dgetrf_(&n,&n,mat,&n,pivot,&fail);
return fail;
}
int nleq_linsol(int n, double *b, struct NLEQ_OPT *opt)
{ int fail=0, nrhs=1;
char mode=JMODE;
if (mat) dgetrs_(&mode,&n,&nrhs,mat,&n,pivot,b,&n,&fail);
else fail = -980;
return fail;
}
void nleq_linfree()
{ if(pivot) free(pivot); pivot = NULL;
if (v) free(v); v = NULL;
if (mat) free(mat); mat = NULL;
}
int nleq_numjac(NLEQ_FFUN *f, int n, double *x, double *fx,
double *scale, JACOBIAN *jac, int *nfcn, struct NLEQ_OPT *opt)
{ int i, k, kn;
int fail=0;
double w, u;
if (!v) fail=nleq_fwalloc(n,&v,"v");
if (fail!=0) return -992;
for (k=0;k