|
Packit |
67cb25 |
/* ode-initval2/modnewton1.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2008, 2009, 2010 Tuomo Keskitalo
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* A modified Newton iteration method for solving non-linear
|
|
Packit |
67cb25 |
equations in implicit Runge-Kutta methods.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* References:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Ascher, U.M., Petzold, L.R., Computer methods for ordinary
|
|
Packit |
67cb25 |
differential and differential-algebraic equations, SIAM,
|
|
Packit |
67cb25 |
Philadelphia, 1998. ISBN 0898714125, 9780898714128
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Hairer, E., Wanner, G., Solving Ordinary Differential
|
|
Packit |
67cb25 |
Equations II: Stiff and Differential-Algebraic Problems,
|
|
Packit |
67cb25 |
Springer, 1996. ISBN 3540604529, 9783540604525
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <string.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_odeiv2.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "odeiv_util.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
typedef struct
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* iteration matrix I - h A (*) J */
|
|
Packit |
67cb25 |
gsl_matrix *IhAJ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* permutation for LU-decomposition */
|
|
Packit |
67cb25 |
gsl_permutation *p;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* difference vector for kth Newton iteration */
|
|
Packit |
67cb25 |
gsl_vector *dYk;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* scaled dYk with desired error level */
|
|
Packit |
67cb25 |
gsl_vector *dScal;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* current Runge-Kutta points (absolute values) */
|
|
Packit |
67cb25 |
double *Yk;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* function values at points Yk */
|
|
Packit |
67cb25 |
double *fYk;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* right hand side of Newton iteration formula */
|
|
Packit |
67cb25 |
gsl_vector *rhs;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* stopping criterion value from previous step */
|
|
Packit |
67cb25 |
double eeta_prev;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
modnewton1_state_t;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void *
|
|
Packit |
67cb25 |
modnewton1_alloc (size_t dim, size_t stage)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
modnewton1_state_t *state =
|
|
Packit |
67cb25 |
(modnewton1_state_t *) malloc (sizeof (modnewton1_state_t));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for modnewton1_state_t",
|
|
Packit |
67cb25 |
GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->IhAJ = gsl_matrix_alloc (dim * stage, dim * stage);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->IhAJ == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (state);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for IhAJ", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->p = gsl_permutation_alloc (dim * stage);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->p == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free (state->IhAJ);
|
|
Packit |
67cb25 |
free (state);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for p", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->dYk = gsl_vector_alloc (dim * stage);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->dYk == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_permutation_free (state->p);
|
|
Packit |
67cb25 |
gsl_matrix_free (state->IhAJ);
|
|
Packit |
67cb25 |
free (state);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for dYk", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->dScal = gsl_vector_alloc (dim * stage);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->dScal == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_free (state->dYk);
|
|
Packit |
67cb25 |
gsl_permutation_free (state->p);
|
|
Packit |
67cb25 |
gsl_matrix_free (state->IhAJ);
|
|
Packit |
67cb25 |
free (state);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for dScal", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->Yk = (double *) malloc (dim * stage * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->Yk == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_free (state->dScal);
|
|
Packit |
67cb25 |
gsl_vector_free (state->dYk);
|
|
Packit |
67cb25 |
gsl_permutation_free (state->p);
|
|
Packit |
67cb25 |
gsl_matrix_free (state->IhAJ);
|
|
Packit |
67cb25 |
free (state);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for Yk", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->fYk = (double *) malloc (dim * stage * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->fYk == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (state->Yk);
|
|
Packit |
67cb25 |
gsl_vector_free (state->dScal);
|
|
Packit |
67cb25 |
gsl_vector_free (state->dYk);
|
|
Packit |
67cb25 |
gsl_permutation_free (state->p);
|
|
Packit |
67cb25 |
gsl_matrix_free (state->IhAJ);
|
|
Packit |
67cb25 |
free (state);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for Yk", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->rhs = gsl_vector_alloc (dim * stage);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->rhs == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (state->fYk);
|
|
Packit |
67cb25 |
free (state->Yk);
|
|
Packit |
67cb25 |
gsl_vector_free (state->dScal);
|
|
Packit |
67cb25 |
gsl_vector_free (state->dYk);
|
|
Packit |
67cb25 |
gsl_permutation_free (state->p);
|
|
Packit |
67cb25 |
gsl_matrix_free (state->IhAJ);
|
|
Packit |
67cb25 |
free (state);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("failed to allocate space for rhs", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->eeta_prev = GSL_DBL_MAX;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return state;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
modnewton1_init (void *vstate, const gsl_matrix * A,
|
|
Packit |
67cb25 |
const double h, const gsl_matrix * dfdy,
|
|
Packit |
67cb25 |
const gsl_odeiv2_system * sys)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Initializes the method by forming the iteration matrix IhAJ
|
|
Packit |
67cb25 |
and generating its LU-decomposition
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
modnewton1_state_t *state = (modnewton1_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix *const IhAJ = state->IhAJ;
|
|
Packit |
67cb25 |
gsl_permutation *const p = state->p;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const size_t dim = sys->dimension;
|
|
Packit |
67cb25 |
const size_t stage = A->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->eeta_prev = GSL_DBL_MAX;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Generate IhAJ */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j, k, m;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < dim; i++)
|
|
Packit |
67cb25 |
for (j = 0; j < dim; j++)
|
|
Packit |
67cb25 |
for (k = 0; k < stage; k++)
|
|
Packit |
67cb25 |
for (m = 0; m < stage; m++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t x = dim * k + i;
|
|
Packit |
67cb25 |
size_t y = dim * m + j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (x != y)
|
|
Packit |
67cb25 |
gsl_matrix_set (IhAJ, x, y,
|
|
Packit |
67cb25 |
-h * gsl_matrix_get (A, k, m) *
|
|
Packit |
67cb25 |
gsl_matrix_get (dfdy, i, j));
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
gsl_matrix_set (IhAJ, x, y,
|
|
Packit |
67cb25 |
1.0 - h * gsl_matrix_get (A, k, m) *
|
|
Packit |
67cb25 |
gsl_matrix_get (dfdy, i, j));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* decompose */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int signum;
|
|
Packit |
67cb25 |
int s = gsl_linalg_LU_decomp (IhAJ, p, &signum);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s != GSL_SUCCESS)
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
modnewton1_solve (void *vstate, const gsl_matrix * A,
|
|
Packit |
67cb25 |
const double c[], const double t, const double h,
|
|
Packit |
67cb25 |
const double y0[], const gsl_odeiv2_system * sys,
|
|
Packit |
67cb25 |
double YZ[], const double errlev[])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Solves the non-linear equation system resulting from implicit
|
|
Packit |
67cb25 |
Runge-Kutta methods by a modified Newton iteration. The
|
|
Packit |
67cb25 |
modified Newton iteration with this problem is stated in the
|
|
Packit |
67cb25 |
form
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
IhAJ * dYk = rhs
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
in which IhAJ is the iteration matrix: IhAJ = (I - hA (*) J) in
|
|
Packit |
67cb25 |
which (*) is the Kronecker matrix product (size of IhAJ =
|
|
Packit |
67cb25 |
dim*stage, dim*stage), dYk is the Runge-Kutta point (Y)
|
|
Packit |
67cb25 |
difference vector for kth Newton iteration: dYk = Y(k+1) - Y(k),
|
|
Packit |
67cb25 |
and rhs = Y(k) - y0 - h * sum j=1..stage (a_j * f(Y(k)))
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves dYk by LU-decomposition of IhAJ with partial
|
|
Packit |
67cb25 |
pivoting.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
modnewton1_state_t *state = (modnewton1_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix *const IhAJ = state->IhAJ;
|
|
Packit |
67cb25 |
gsl_permutation *const p = state->p;
|
|
Packit |
67cb25 |
gsl_vector *const dYk = state->dYk;
|
|
Packit |
67cb25 |
double *const Yk = state->Yk;
|
|
Packit |
67cb25 |
double *const fYk = state->fYk;
|
|
Packit |
67cb25 |
gsl_vector *const rhs = state->rhs;
|
|
Packit |
67cb25 |
double *const eeta_prev = &(state->eeta_prev);
|
|
Packit |
67cb25 |
gsl_vector *const dScal = state->dScal;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const size_t dim = sys->dimension;
|
|
Packit |
67cb25 |
const size_t stage = A->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int status = GSL_CONTINUE; /* Convergence status for Newton iteration */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Set starting values for iteration. Use starting values of Y(k) =
|
|
Packit |
67cb25 |
y0. FIXME: Implement better initial guess described in Hairer and
|
|
Packit |
67cb25 |
Wanner.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t j, m;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set_zero (dYk);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < stage; j++)
|
|
Packit |
67cb25 |
for (m = 0; m < dim; m++)
|
|
Packit |
67cb25 |
Yk[j * dim + m] = y0[m];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Loop modified Newton iterations */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t iter = 0;
|
|
Packit |
67cb25 |
size_t j, m, n;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Maximum number of Newton iterations. */
|
|
Packit |
67cb25 |
const size_t max_iter = 7;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double dScal_norm = 0.0; /* Newton iteration step length */
|
|
Packit |
67cb25 |
double dScal_norm_prev = 0.0; /* Previous dScal_norm */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while (status == GSL_CONTINUE && iter < max_iter)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
iter++;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Update Y(k) and evaluate function */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < stage; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (m = 0; m < dim; m++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
Yk[j * dim + m] += gsl_vector_get (dYk, j * dim + m);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = GSL_ODEIV_FN_EVAL (sys, t + c[j] * h, &Yk[j * dim],
|
|
Packit |
67cb25 |
&fYk[j * dim]);
|
|
Packit |
67cb25 |
if (s != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Calculate rhs */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < stage; j++)
|
|
Packit |
67cb25 |
for (m = 0; m < dim; m++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double sum = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (n = 0; n < stage; n++)
|
|
Packit |
67cb25 |
sum += gsl_matrix_get (A, j, n) * fYk[n * dim + m];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set (rhs, j * dim + m,
|
|
Packit |
67cb25 |
-1.0 * Yk[j * dim + m] + y0[m] + h * sum);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solve dYk */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = gsl_linalg_LU_solve (IhAJ, p, rhs, dYk);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Evaluate convergence according to method by Hairer and
|
|
Packit |
67cb25 |
Wanner, section IV.8. The iteration step is normalized by
|
|
Packit |
67cb25 |
desired error level before calculation of the norm, to take
|
|
Packit |
67cb25 |
the tolerance on each component of y into account.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* relative change in two Newton iteration steps */
|
|
Packit |
67cb25 |
double theta_k;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* transformation of theta_k */
|
|
Packit |
67cb25 |
double eeta_k = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < stage; j++)
|
|
Packit |
67cb25 |
for (m = 0; m < dim; m++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_set (dScal, j * dim + m,
|
|
Packit |
67cb25 |
gsl_vector_get (dYk, j * dim + m)
|
|
Packit |
67cb25 |
/ errlev[m]);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
dScal_norm_prev = dScal_norm;
|
|
Packit |
67cb25 |
dScal_norm = gsl_blas_dnrm2 (dScal);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
theta_k = dScal_norm / dScal_norm_prev;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (iter > 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* check for increase in step size, which indicates divergence */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (theta_k >= 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_FAILURE;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
eeta_k = theta_k / (1.0 - theta_k);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
eeta_k = pow (GSL_MAX_DBL (*eeta_prev, GSL_DBL_EPSILON), 0.8);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*eeta_prev = eeta_k;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (eeta_k * dScal_norm < 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
status = GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* No convergence reached within allowed iterations */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_FAILURE;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Give solution in YZ */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t j, m;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < stage; j++)
|
|
Packit |
67cb25 |
for (m = 0; m < dim; m++)
|
|
Packit |
67cb25 |
YZ[j * dim + m] =
|
|
Packit |
67cb25 |
Yk[j * dim + m] + gsl_vector_get (dYk, j * dim + m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
modnewton1_free (void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
modnewton1_state_t *state = (modnewton1_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_free (state->rhs);
|
|
Packit |
67cb25 |
free (state->fYk);
|
|
Packit |
67cb25 |
free (state->Yk);
|
|
Packit |
67cb25 |
gsl_vector_free (state->dScal);
|
|
Packit |
67cb25 |
gsl_vector_free (state->dYk);
|
|
Packit |
67cb25 |
gsl_permutation_free (state->p);
|
|
Packit |
67cb25 |
gsl_matrix_free (state->IhAJ);
|
|
Packit |
67cb25 |
free (state);
|
|
Packit |
67cb25 |
}
|