# «NITSOL: A NEWTON ITERATIVE SOLVER FOR NONLINEAR SYSTEMS∗ MICHAEL PERNICE† AND HOMER F. WALKER‡ Abstract. We introduce a well-developed Newton ...»

c 1998 Society for Industrial and Applied Mathematics

SIAM J. SCI. COMPUT.

Vol. 19, No. 1, pp. 302–318, January 1998 021

## NITSOL: A NEWTON ITERATIVE SOLVER

## FOR NONLINEAR SYSTEMS∗

## MICHAEL PERNICE† AND HOMER F. WALKER‡

Abstract. We introduce a well-developed Newton iterative (truncated Newton) algorithm for solving large-scale nonlinear systems. The framework is an inexact Newton method globalized by backtracking. Trial steps are obtained using one of several Krylov subspace methods. The algorithm is implemented in a Fortran solver called NITSOL that is robust yet easy to use and provides a number of useful options and features. The structure oﬀers the user great ﬂexibility in addressing problem speciﬁcity through preconditioning and other means and allows easy adaptation to parallel environments. Features and capabilities are illustrated in numerical experiments.

Key words. Newton iterative methods, truncated Newton methods, Newton–Krylov methods, inexact Newton methods, Newton’s method, forcing terms, iterative linear algebra methods, Krylov subspace methods, GMRES, BiCGSTAB, TFQMR AMS subject classiﬁcations. 65H10, 65F10, 65N22 PII. S1064827596303843

1. Introduction. Our interest is in methods for solving a system of nonlinear equations F : Rn → Rn, (1.1) F (x) = 0, where F is assumed to be continuously diﬀerentiable everywhere in Rn. Newton’s method for solving (1.1) requires, at the kth step, the solution of the linear Newton equation F (xk ) sk = −F (xk ), (1.2) where xk is the current approximate solution. A Newton iterative method, or truncated Newton method, is an implementation of Newton’s method in which an iterative linear solver is used to determine an approximate solution of (1.2). Newton iterative methods are especially well suited for large-scale problems and have been used very successfully in many scientiﬁc and industrial applications.

Our purpose here is to introduce a well-developed Newton iterative algorithm and to describe its implementation in a Fortran solver called NITSOL. In brief, the underlying nonlinear algorithm is an inexact Newton method with a backtracking globalization1, previously formulated as Algorithm INB of [9, section 6]. In this, the Newton equation (1.2) is relaxed to an inexact Newton condition (cf. [8]) F (xk ) + F (xk ) sk ≤ ηk F (xk ), (1.3) ∗ Received by the editors May 17, 1996; accepted for publication (in revised form) March 10, 1997.

http://www.siam.org/journals/sisc/19-1/30384.html † Center for High Performance Computing, University of Utah, Salt Lake City, UT 84112 (usimap@sneﬀels.usi.utah.edu).

‡ Department of Mathematics and Statistics, Utah State University, Logan, UT 84322-3900 (walker@math.usu.edu). Current address: Mathematical Sciences Department, Worcester Polytechnic Institute, Worcester, MA, 01609. The research of this author was supported in part by U.S.

Department of Energy grant DE-FG03-94ER25221 and National Science Foundation grant DMSboth with Utah State University.

1 We use the term “globalization” in its traditional sense, i.e., to refer to strategies for improving

in which the “forcing term” ηk ∈ [0, 1) can be speciﬁed in several ways as in [10] to enhance eﬃciency and convergence. An initial sk satisfying (1.3) is determined by using a Krylov subspace method to solve (1.2) approximately; thus, our algorithm is a Newton–Krylov method. The Krylov subspace methods implemented in NITSOL are GMRES(m) [25], [33], BiCGSTAB [31], and TFQMR [12]. Once an initial sk has been determined, it is tested and, if necessary, reduced in length through safeguarded backtracking until an acceptable step is obtained.

With NITSOL, we have hoped to provide a robust, theoretically well-founded solver that is simple and easy to use but which oﬀers enough options and ﬂexibility to allow the user to implement sophisticated solution strategies that address speciﬁc problem needs and particular machine capabilities. No preconditioners are provided by NITSOL, but the structure allows the user to implement a preferred preconditioner with minimal diﬃculty. The user is also allowed to specify or supply an inner product and associated norm that are used throughout the algorithm; this capability not only provides a means of addressing problem scaling but also allows easy adaptation to parallel environments. There are options for producing a variety of diagnostic information and also for passing to user-supplied routines information that can be helpful in determining when to re-evaluate Jacobians or preconditioners.

Several Newton iterative solvers have been introduced previously by others. The NKSOL Fortran code [3] implements a Newton iterative method with either backtracking or trust region globalization; the Krylov solver options are GMRES and the full orthogonalization (Arnoldi) method [23]. The Scalable Nonlinear Equations Solver (SNES) package [15], which is part of the very extensive C library PETSc [16], oﬀers a Newton iterative solver with either backtracking or trust region globalization and provides various Krylov solvers and preconditioners. SNES and PETSc achieve parallelism through message passing, using the Message-Passing Interface standard (MPI) [28] for interprocessor communication. The more specialized parallel reactive ﬂow code MPSalsa [26], [27], also written in C, includes a nonlinear solver that, like NITSOL, is based on Algorithm INB of [9]; various Krylov solvers and incomplete factorization preconditioners are provided via the Aztec parallel iterative linear algebra library [17]. Less closely related are the truncated Newton solvers LANCELOT [6] and TN [21], [22], which employ the conjugate gradient method and are intended primarily for optimization.

In section 2 below, we introduce the inexact Newton backtracking method, discuss its implementation as a Newton iterative method, and outline the Krylov solver options. In section 3, we discuss the structure, features, and usage of NITSOL. In section 4, we give illustrative examples of NITSOL applications and usage.

2. The algorithm. The inexact Newton backtracking method that we implement is the following from [9], which oﬀers strong global convergence properties combined with potentially fast local convergence. In this, the inexact Newton condition (1.3) is augmented with a condition of suﬃcient reduction of the norm of F.

ALGORITHM INB (Inexact Newton Backtracking method [9]).

Let x0, ηmax ∈ [0, 1), t ∈ (0, 1), and 0 θmin θmax 1 be given.

**For k = 0, 1,... (until convergence) do:**

Choose initial ηk ∈ [0, ηmax ] and sk such that F (xk ) + F (xk ) sk ≤ ηk F (xk ).

**While F (xk + sk ) [1 − t(1 − ηk )] F (xk ) do:**

Choose θ ∈ [θmin, θmax ].

304 MICHAEL PERNICE AND HOMER F. WALKER Update sk ←− θsk and ηk ←− 1 − θ(1 − ηk ).

Set xk+1 = xk + sk.

Note that, given an initial ηk, a satisfactory initial sk exists if the Newton equation (1.2) is consistent, in particular, if F (xk ) is invertible. If a satisfactory initial sk can be found, then we have from remarks in [9, section 6] that in exact arithmetic Algorithm INB does not break down in the while-loop, i.e., an acceptable sk is determined after, at most, a ﬁnite number of step reductions. Furthermore, it is easy to see that an inexact Newton condition (1.3) holds for each sk and ηk determined by the while-loop and, in particular, for the ﬁnal sk and ηk.

The principal theoretical result for Algorithm INB in exact arithmetic is the following.

** THEOREM 2.1 (see [9]).**

Assume that F is continuously diﬀerentiable. If {xk } produced by Algorithm INB has a limit point x∗ such that F (x∗ ) is invertible, then F (x∗ ) = 0 and xk → x∗. Furthermore, the initial sk and ηk are accepted without modiﬁcation in the while-loop for all suﬃciently large k.

It follows that if Algorithm INB does not terminate and produces {xk }, then

**exactly one of the following must hold:**

• xk → ∞, i.e., {xk } has no limit points.

• {xk } has one or more limit points, and F is singular at each of them.

• {xk } converges to a solution x∗ at which F is invertible.

Note that in the case of the (desirable) third alternative, the initial sk and ηk are accepted without modiﬁcation for all suﬃciently large k and, consequently, asymptotic convergence to the solution is determined by the initial ηk ’s as in the local convergence analysis of [8].

To implement Algorithm INB as a Newton iterative method, one must ﬁrst specify various details in the algorithm and then choose a suitable iterative linear solver for determining initial inexact Newton steps. In our implementation, relatively minor details are speciﬁed as follows: as in [10, section 3.1], we use backtracking safeguard values θmin = 0.1 and θmax = 0.5. We take the norm to be an inner-product norm and choose θ ∈ [θmin, θmax ] to minimize a quadratic that interpolates F in the direction of the inexact Newton step. The value t = 10−4 is used to judge suﬃcient reduction, and we use the upper bound ηmax = 0.9. Convergence is declared if either F (xk ) ≤ FTOL or sk ≤ STPTOL · xk, where FTOL and STPTOL are user-supplied tolerances.

In the remainder of this section, we ﬁrst address a major detail: choosing the initial forcing terms. We then discuss the Krylov solvers that are available in our implementation.

2.1. Choosing the forcing terms. As noted above, if the iterates produced by Algorithm INB converge to a solution of (1.1) at which F is invertible, then the ultimate speed of convergence is determined by the initial forcing terms ηk. In particular, by the analysis of [8], desirably fast local convergence can be obtained by making suitably small choices of the initial forcing terms near a solution.

The initial forcing terms can also signiﬁcantly aﬀect the performance of the algorithm away from a solution. Indeed, choosing an initial ηk too small can lead to oversolving the Newton equation, i.e., imposing an accuracy on an approximate solution sk of (1.2) that leads to signiﬁcant disagreement between F (xk + sk ) and its local linear model F (xk ) + F (xk )sk. Oversolving may result in little or no progress toward a solution. Moreover, in a Newton iterative implementation, it may involve pointless expense; a less accurate solution of (1.2) may be both cheaper and more eﬀective in reducing the norm of F.

305

## NITSOL: A NEWTON ITERATIVE SOLVER

Here, we implement two choices of the initial forcing terms from [10] that give desirably fast local convergence and also tend to minimize oversolving. These are as**follows:**

Choice 1: Select any η0 ∈ [0, 1) and choose

Also, in our implementation, we use the initial value η0 = 0.5 with both Choice 1 and Choice 2.

Theorems 2.2 and 2.3 in [10] give results concerning the speed of local convergence to a solution of (1.1) when (2.1) or (2.2) is used to determine inexact Newton iterates.

Augmenting Theorem 2.1 above with those results immediately yields the following extensions.

** THEOREM 2.2.**

Assume that F is continuously diﬀerentiable and that each ηk in Algorithm INB is given by (2.1) followed by (2.3). If {xk } produced by Algorithm INB has a limit point x∗ such that F (x∗ ) is invertible, then F (x∗ ) = 0 and xk → x∗.

Furthermore, if F is Lipschitz continuous at x∗, then

We employ one further safeguard with the above choices that plays a role near a solution of (1.1) when a stopping criterion of the form F (xk ) is used. (Here, represents either an absolute or a relative tolerance.) Maintaining agreement between F and its local linear model means that

** F (xk+1 ) ≈ F (xk ) + F (xk )sk ≤ ηk F (xk ).**

Thus, if ηk F (xk ), then it is likely that the stopping criterion can be satisﬁed with less eﬀort by using a larger ηk. On the other hand, if ηk F (xk ) is slightly larger than, then it is possible that the stopping criterion can be satisﬁed without an additional inexact Newton step by using a slightly smaller ηk. Therefore, we impose the following ﬁnal safeguard: if ηk ≤ 2 / F (xk ), then ηk ←− 0.8 / F (xk ).

In addition to the above choices, we also allow the user to specify a constant choice ηk = η∗, where η∗ ∈ [0, 1) is independent of k. Such a choice is occasionally useful, e.g., when only modest accuracy is desired in approximately solving (1.2). In this case, augmenting Theorem 2.1 above with Theorem 2.3 in [8] immediately gives the following convergence result for this choice.

** THEOREM 2.4.**

Assume that F is continuously diﬀerentiable and that each ηk in Algorithm INB is given by ηk = η∗ for some η∗ ∈ [0, 1) independent of k. If {xk } produced by Algorithm INB has a limit point x∗ such that F (x∗ ) is invertible, then F (x∗ ) = 0 and xk → x∗. Furthermore, for every η ∈ (η∗, 1), we have

2.2. The Krylov subspace methods. In our Newton iterative implementation of Algorithm INB, a Krylov subspace method is used to obtain each initial sk. Krylov subspace methods comprise a very broad and successful class of iterative linear algebra methods; see [13] for a recent survey.

As indicated in section 1, the Krylov solver options in our algorithm are GMRES(m), BiCGSTAB, and TFQMR. These methods are well known, widely used, and applicable to general linear systems. They are “transpose free”, i.e., in an implementation, they require only products of the coeﬃcient matrix with vectors and do not require products involving its transpose. Being “transpose free” can be very advantageous, particularly in the Newton iterative context, in which the coeﬃcient matrix is F (xk ). If analytic evaluation of products of F (xk ) with vectors is undesirable or infeasible, then these products can be approximated with ﬁnite diﬀerences of F -values, as discussed further below.