From f6a2235bd198c3c53746ff698e3ca52d3ae09577 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Tue, 17 Jan 2023 22:38:40 -0600 Subject: [PATCH] added files from https://github.com/NCAR/ncl added readme and VSCode workspace file --- LICENSE | 13 ++ README.md | 5 + splpak.code-workspace | 13 ++ src/bascmpd.f | 221 +++++++++++++++++++++ src/cfaerr.f | 40 ++++ src/splccd.f | 406 ++++++++++++++++++++++++++++++++++++++ src/splcwd.f | 439 ++++++++++++++++++++++++++++++++++++++++++ src/splded.f | 116 +++++++++++ src/splfed.f | 27 +++ src/suprld.f | 288 +++++++++++++++++++++++++++ 10 files changed, 1568 insertions(+) create mode 100644 LICENSE create mode 100644 README.md create mode 100644 splpak.code-workspace create mode 100644 src/bascmpd.f create mode 100644 src/cfaerr.f create mode 100644 src/splccd.f create mode 100644 src/splcwd.f create mode 100644 src/splded.f create mode 100644 src/splfed.f create mode 100644 src/suprld.f diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..f2f5211 --- /dev/null +++ b/LICENSE @@ -0,0 +1,13 @@ +Copyright 2019 University Corporation for Atmospheric Research + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. diff --git a/README.md b/README.md new file mode 100644 index 0000000..ef16891 --- /dev/null +++ b/README.md @@ -0,0 +1,5 @@ +Double precision SPLPAK files from [NCL](https://github.com/NCAR/ncl). + +### See also + * [bspline-fortran](https://github.com/jacobwilliams/bspline-fortran) Multidimensional B-Spline Interpolation of Data on a Regular Grid + * [bspline](https://github.com/NCAR/bspline) - Cubic B-Spline implementation in C++ templates. Also has a copy of [splpak.f](https://github.com/NCAR/bspline/tree/master/Tests/Fortran) diff --git a/splpak.code-workspace b/splpak.code-workspace new file mode 100644 index 0000000..489f23c --- /dev/null +++ b/splpak.code-workspace @@ -0,0 +1,13 @@ +{ + "folders": [ + { + "path": "." + } + ], + "settings": { + "files.trimTrailingWhitespace": true, + "editor.insertSpaces": true, + "editor.tabSize": 4, + "editor.trimAutoWhitespace": true + } +} \ No newline at end of file diff --git a/src/bascmpd.f b/src/bascmpd.f new file mode 100644 index 0000000..07f099a --- /dev/null +++ b/src/bascmpd.f @@ -0,0 +1,221 @@ +C +C Copyright (C) 2000 +C University Corporation for Atmospheric Research +C All Rights Reserved +C +C The use of this Software is governed by a License Agreement. +C + SUBROUTINE BASCMPD(X,NDERIV,XMIN,NODES,ICOL,BASM) + DOUBLE PRECISION X + DOUBLE PRECISION XMIN + DOUBLE PRECISION BASM + DOUBLE PRECISION DX + DOUBLE PRECISION DXIN + DOUBLE PRECISION XB + DOUBLE PRECISION BAS1 + DOUBLE PRECISION Z + DOUBLE PRECISION FACT + DOUBLE PRECISION Z1 +C +C This routine does basis function computations for natural +C splines. This routine is called by routines SPLCW and SPLDE +C to compute ICOL and BASM, which are defined as follows: +C +C The MDIM indices in IB (defined through common) determine +C a specific node in the node grid (see routine SPLCC for a +C description of the node grid). Every node is associated +C with an MDIM-dimensional basis function and a corresponding +C column in the least squares matrix (or element of the +C coefficient vector). The column index (which may be thought +C of as a linear address for the MDIM-dimensional node grid) +C corresponding to the specified node is computed as ICOL. The +C associated basis function evaluated at X (an MDIM-vector) is +C computed as BASM (a scalar). +C +C In case NDERIV is not all zero, BASM will not be the value of +C the basis function but rather a partial derivative of that +C function as follows: +C +C The order of the partial derivative in the direction of the +C IDIM coordinate is NDERIV(IDIM) (for IDIM .LE. MDIM). This +C routine will compute incorrect values if NDERIV(IDIM) is not +C in the range 0 to 2. +C +C + DIMENSION X(4),NDERIV(4),XMIN(4),NODES(4) +C +C The technique of this routine is to transform the independent +C variable in each dimension such that the nodes fall on +C suitably chosen integers. On this transformed space, the +C 1-dimensional basis functions and their derivatives have a +C particularly simple form. The desired MDIM-dimensional basis +C function (or any of its partial derivatives) is computed as +C a product of such 1-dimensional functions (tensor product +C method of defining multi-dimensional splines). The values +C which determine the location of the nodes, and hence the +C above transform, are passed through common and the argument +C list. +C + COMMON /SPLCOMD/DX(4),DXIN(4),MDIM,IB(4),IBMN(4),IBMX(4) + SAVE +C +C ICOL will be a linear address corresponding to the indices in IB. +C + ICOL = 0 +C +C BASM will be M-dimensional basis function evaluated at X. +C + BASM = 1.D0 + DO 121 IDIM = 1,MDIM +C +C Compute ICOL by Horner's method. +C + MDMID = MDIM + 1 - IDIM + ICOL = NODES(MDMID)*ICOL + IB(MDMID) +C +C NGO depends upon function type and NDERIV. +C + NTYP = 1 +C +C Function type 1 (left linear) for IB = 0 or 1. +C + IF (IB(IDIM).LE.1) GO TO 101 + NTYP = 2 +C +C Function type 2 (chapeau function) for 2 LT IB LT NODES-2. +C + IF (IB(IDIM).LT.NODES(IDIM)-2) GO TO 101 + NTYP = 3 +C +C Function type 3 (right linear) for IB = NODES-2 or NODES-1. +C + 101 NGO = 3*NTYP + NDERIV(IDIM) - 2 +C +C XB is X value of node IB (center of basis function). +C + XB = XMIN(IDIM) + DBLE(IB(IDIM))*DX(IDIM) +C +C BAS1 will be the 1-dimensional basis function evaluated at X. +C + BAS1 = 0.D0 + GO TO (102,103,104,105,106,108,110,113,117) NGO +C +C Function type 1 (left linear) is mirror image of function type 3. +C +C Transform so that XB is at 2 and the other nodes are at the integers +C (with ordering reversed to form a mirror image). +C + 102 Z = DXIN(IDIM)* (XB-X(IDIM)) + 2.D0 + GO TO 111 +C +C 1st derivative. +C + 103 FACT = -DXIN(IDIM) + GO TO 114 +C +C 2nd derivative. +C + 104 FACT = -DXIN(IDIM) + GO TO 118 +C +C Function type 2 (chapeau function). +C +C Transform so that XB is at the origin and the other nodes are at +C the integers. +C + 105 Z = ABS(DXIN(IDIM)* (X(IDIM)-XB)) - 2.D0 +C +C This chapeau function is then that unique cubic spline which is +C identically zero for ABS(Z) GE 2 and is 1 at the origin. This +C function is the general interior node basis function. +C + IF (Z.GE.0.D0) GO TO 120 + BAS1 = -.25D0*Z**3 + Z = Z + 1.D0 + IF (Z.GE.0.D0) GO TO 120 + BAS1 = BAS1 + Z**3 + GO TO 120 +C +C 1st derivative. +C + 106 Z = X(IDIM) - XB + FACT = DXIN(IDIM) + IF (Z.LT.0.D0) FACT = -FACT + Z = FACT*Z - 2.D0 + IF (Z.GE.0.D0) GO TO 120 + BAS1 = -.75D0*Z**2 + Z = Z + 1.D0 + IF (Z.GE.0.D0) GO TO 107 + BAS1 = BAS1 + 3.D0*Z**2 + 107 BAS1 = FACT*BAS1 + GO TO 120 +C +C 2nd derivative. +C + 108 FACT = DXIN(IDIM) + Z = FACT*ABS(X(IDIM)-XB) - 2.D0 + IF (Z.GE.0.D0) GO TO 120 + BAS1 = -1.5D0*Z + Z = Z + 1.D0 + IF (Z.GE.0.D0) GO TO 109 + BAS1 = BAS1 + 6.D0*Z + 109 BAS1 = (FACT**2)*BAS1 + GO TO 120 +C +C Function type 3 (right linear). +C +C Transform so that XB is at 2 and the other nodes are at the integers. +C + 110 Z = DXIN(IDIM)* (X(IDIM)-XB) + 2.D0 +C +C This right linear function is defined to be that unique cubic spline +C which is identically zero for Z .LE. 0 and is 3*Z-3 for Z GE 2. +C This function (obviously having zero 2nd derivative for +C Z GE 2) is used for the two nodes nearest an edge in order +C to generate natural splines, which must by definition have +C zero 2nd derivative at the boundary. +C +C Note that this method of generating natural splines also provides +C a linear extrapolation which has 2nd order continuity with +C the interior splines at the boundary. +C + 111 IF (Z.LE.0.D0) GO TO 120 + IF (Z.GE.2.D0) GO TO 112 + BAS1 = .5D0*Z**3 + Z = Z - 1.D0 + IF (Z.LE.0.D0) GO TO 120 + BAS1 = BAS1 - Z**3 + GO TO 120 + 112 BAS1 = 3.D0*Z - 3.D0 + GO TO 120 +C +C 1st derivative. +C + 113 FACT = DXIN(IDIM) + 114 Z = FACT* (X(IDIM)-XB) + 2.D0 + IF (Z.LE.0.D0) GO TO 120 + IF (Z.GE.2.D0) GO TO 116 + BAS1 = 1.5D0*Z**2 + Z = Z - 1.D0 + IF (Z.LE.0.D0) GO TO 115 + BAS1 = BAS1 - 3.D0*Z**2 + 115 BAS1 = FACT*BAS1 + GO TO 120 + 116 BAS1 = 3.D0*FACT + GO TO 120 +C +C 2nd derivative. +C + 117 FACT = DXIN(IDIM) + 118 Z = FACT* (X(IDIM)-XB) + 2.D0 + Z1 = Z - 1.D0 + IF (ABS(Z1).GE.1.D0) GO TO 120 + BAS1 = 3.D0*Z + IF (Z1.LE.0.D0) GO TO 119 + BAS1 = BAS1 - 6.D0*Z1 + 119 BAS1 = (FACT**2)*BAS1 + 120 BASM = BASM*BAS1 + 121 CONTINUE + ICOL = ICOL + 1 + RETURN + END diff --git a/src/cfaerr.f b/src/cfaerr.f new file mode 100644 index 0000000..e362b91 --- /dev/null +++ b/src/cfaerr.f @@ -0,0 +1,40 @@ +C +C Copyright (C) 2000 +C University Corporation for Atmospheric Research +C All Rights Reserved +C +C The use of this Software is governed by a License Agreement. +C +C SUBROUTINE CFAERR (IERR,MESS,LMESS) +C +C PURPOSE To print an error number and an error message +C or just an error message. +C +C USAGE CALL CFAERR (IERR,MESS,LMESS) +C +C ARGUMENTS +C ON INPUT IERR +C The error number (printed only if non-zero). +C +C MESS +C Message to be printed. +C +C LMESS +C Number of characters in mess (.LE. 130). +C +C ARGUMENTS +C ON OUTPUT None +C +C I/O The message is writen to unit 6. +C +C ****************************************************************** +C + SUBROUTINE CFAERR (IERR,MESS,LMESS) +C + CHARACTER *(*) MESS +C + IF (IERR .NE. 0) WRITE (6,'(A,I5)') ' IERR=', IERR + WRITE (6,'(A)') MESS(1:LMESS) +C + RETURN + END diff --git a/src/splccd.f b/src/splccd.f new file mode 100644 index 0000000..76b9d56 --- /dev/null +++ b/src/splccd.f @@ -0,0 +1,406 @@ +C PACKAGE SPLPAK Documentation for user entries follows +C the general package information. +C +C LATEST REVISION August, 1998 +C +C PURPOSE This package contains routines for fitting +C (least squares) a multidimensional cubic spline +C to arbitrarily located data. It also contains +C routines for evaluating this spline (or its +C partial derivatives) at any point. +C +C Coefficient calculation is performed in +C subroutines SPLCC or SPLCW and evaluation is +C performed by functions SPLFE or SPLDE. +C +C USAGE Package SPLPAK contains four user entries -- +C SPLCC, SPLCW, SPLFE, AND SPLDE. +C +C The user first calls SPLCC by +C +C CALL SPLCC (NDIM,XDATA,L1XDAT,YDATA,NDATA, +C XMIN,XMAX,NODES,XTRAP,COEF,NCF, +C WORK,NWRK,IERROR) +C +C or SPLCW by +C +C CALL SPLCW (NDIM,XDATA,L1XDATA,YDATA,WDATA, +C NDATA,XMIN,XMAX,NODES,XTRAP, +C COEF,NCF,WORK,NWRK,IERROR) +C +C The parameter NDATA in the call to SPLCW +C enables the user to weight some of the data +C points more heavily than others. Both +C routines return a set of coefficients in the +C array COEF. These coefficients are +C subsequently used in the computation of +C function values and partial derivatives. +C To compute values on the spline approximation +C the user then calls SPLFE or SPLDE any +C number of times in any order provided that +C the values of the inputs, NDIM, COEF, XMIN, +C XMAX, and NODES, are preserved between calls. +C +C*PL*ERROR* Comment line too long +C SPLFE and SPLDE are called in the following way: +C +C F = SPLFE (NDIM,X,COEF,XMIN,XMAX,NODES, +C IERROR) +C +C or +C +C F = SPLDE (NDIM,X,NDERIV,COEF,XMIN,XMAX, +C NODES,IERROR) +C +C The routine SPLFE returns an interpolated +C value at the point defined by the array X. +C SPLDE affords the user the additional +C capability of calculating an interpolated +C value for one of several partial derivatives +C specified by the array NDERIV. +C +C I/O None, except for error messages printed by +C calls to CFSARR, if an error is detected. +C +C PRECISION Single +C +C REQUIRED LIBRARY SUPRLS, CFAERR +C FILES +C +C LANGUAGE FORTRAN +C +C HISTORY Developed in 1972-73 by NCAR's +C Scientific Computing Division. +C +C Cleaned up and added to the Ngmath library in +C 1998. +C +C PORTABILITY FORTRAN 77 +C +C*********************************************************************** +C +C SUBROUTINE SPLCCD(NDIM,XDATA,L1XDAT,YDATA,NDATA,XMIN,XMAX,NODES, +C XTRAP,COEF,NCF,WORK,NWRK,IERROR) +C +C DIMENSION OF XDATA(NDATA,L1XDAT),YDATA(NDATA),XMIN(NDIM), +C ARGUMENTS XMAX(NDIM),NODES(NDIM),COEF(NCF),WORK(NWRK) +C +C PURPOSE N-dimensional cubic spline coefficient +C calculation by least squares. +C +C USAGE The usage and arguments of this routine are +C identical to those for SPLCW except for the +C omission of the array of weights, WDATA. See +C entry SPLCW description immediately below for a +C complete description. +C +C CALL SPLCCD(NDIM,XDATA,L1XDAT,YDATA,NDATA,XMIN, +C XMAX,NODES,XTRAP,COEF,NCF,WORK, +C NWRK,IERROR) +C +C*********************************************************************** +C +C SUBROUTINE SPLCWD(NDIM,XDATA,L1XDAT,YDATA,WDATA,NDATA,XMIN,XMAX, +C NODES,XTRAP,COEF,NCF,WORK,NWRK,IERROR) +C +C +C DIMENSION OF XDATA(L1XDAT,NDATA),YDATA(NDATA),WDATA(NDATA), +C ARGUMENTS XMIN(NDIM),XMAX(NDIM),NODES(NDIM),COEF(NCF), +C WORK(NWRK) +C +C PURPOSE N-dimensional cubic spline coefficient +C calculation by weighted least squares on +C arbitrarily located data. +C +C A grid of evenly spaced nodes in NDIM space is +C defined by the arguments XMIN, XMAX and NODES. +C A linear basis for the class of natural splines +C on these nodes is formed, and a set of +C corresponding coefficients is computed in the +C array COEF. These coefficients are chosen to +C minimize the weighted sum of squared errors +C between the spline and the arbitrarily located +C data values described by the arguments XDATA, +C YDATA and NDATA. The smoothness of the spline +C in data sparse areas is controlled by the +C argument XTRAP. +C +C NOTE In order to understand the arguments of this +C routine, one should realize that the node grid +C need not bear any particular relation to the +C data points. In the theory of exact-fit +C interpolatory splines, the nodes would in fact +C be data locations, but in this case they serve +C only to define the class of splines from which +C the approximating function is chosen. This +C node grid is a rectangular arrangement of +C points in NDIM space, with the restriction that +C along any coordinate direction the nodes are +C equally spaced. The class of natural splines +C on this grid of nodes (NDIM-cubic splines whose +C 2nd derivatives normal to the boundaries are 0) +C has as many degrees of freedom as the grid has +C nodes. Thus the smoothness or flexibility of +C the splines is determined by the choice of the +C node grid. +C +C USAGE CALL SPLCW (NDIM,XDATA,L1XDAT,YDATA,WDATA, +C NDATA,XMIN,XMAX,NODES,XTRAP,COEF, +C NCF,WORK,NWRK,IERROR) +C +C The spline (or its derivatives) may then be +C evaluated by using function SPLFE (or SPLDE). +C +C ARGUMENTS +C +C ON INPUT NDIM +C The dimensionality of the problem. The +C spline is a function of NDIM variables or +C coordinates and thus a point in the +C independent variable space is an NDIM vector. +C NDIM must be in the range 1 .LE. NDIM .LE. 4. +C +C XDATA +C A collection of locations for the data +C values, i.e., points from the independent +C variable space. This collection is a +C 2-dimensional array whose 1st dimension +C indexes the NDIM coordinates of a given point +C and whose 2nd dimension labels the data +C point. For example, the data point with +C label IDATA is located at the point +C (XDATA(1,IDATA),...,XDATA(NDIM,IDATA)) where +C the elements of this vector are the values of +C the NDIM coordinates. The location, number +C and ordering of the data points is arbitrary. +C The dimension of XDATA is assumed to be +C XDATA(L1XDAT,NDATA). +C +C L1XDAT +C The length of the 1st dimension of XDATA in +C the calling program. L1XDAT must be .GE. +C NDIM. +C +C NOTE: For 1-dimensional problems L1XDAT +C is usually 1. +C +C YDATA +C A collection of data values corresponding to +C the points in XDATA. YDATA(IDATA) is the +C data value associated with the point +C (XDATA(1,IDATA),...,XDATA(NDIM,IDATA)) in the +C independent variable space. The spline whose +C coefficients are computed by this routine +C approximates these data values in the least +C*PL*ERROR* Comment line too long +C squares sense. The dimension is assumed to be +C YDATA(NDATA). +C +C WDATA +C A collection of weights. WDATA(IDATA) is a +C weight associated with the data point +C labelled IDATA. It should be non-negative, +C but may be of any magnitude. The weights +C have the effect of forcing greater or lesser +C accuracy at a given point as follows: this +C routine chooses coefficients to minimize the +C sum over all data points of the quantity +C +C (WDATA(IDATA)*(YDATA(IDATA) - spline value +C at XDATA(IDATA)))**2. +C +C Thus, if the reliability +C of a data point is known to be low, the +C corresponding weight may be made small +C (relative to the other weights) so that the +C sum over all data points is affected less by +C discrepencies at the unreliable point. Data +C points with zero weight are completely +C ignored. +C +C NOTE: If WDATA(1) is .LT. 0, the other +C elements of WDATA are not +C referenced, and all weights are +C assumed to be unity. +C +C The dimension is assumed to be WDATA(NDATA) +C unless WDATA(1) .LT. 0., in which case the +C dimension is assumed to be 1. +C +C NDATA +C The number of data points mentioned in the +C above arguments. +C +C XMIN +C A vector describing the lower extreme corner +C of the node grid. A set of evenly spaced +C nodes is formed along each coordinate axis +C and XMIN(IDIM) is the location of the first +C node along the IDIM axis. The dimension is +C assumed to be XMIN(NDIM). +C +C XMAX +C A vector describing the upper extreme corner +C of the node grid. A set of evenly spaced +C nodes is formed along each coordinate axis +C and XMAX(IDIM) is the location of the last +C node along the IDIM axis. The dimension is +C assumed to be XMAX(NDIM). +C +C NODES +C A vector of integers describing the number of +C nodes along each axis. NODES(IDIM) is the +C number of nodes (counting endpoints) along +C the IDIM axis and determines the flexibility +C of the spline in that coordinate direction. +C NODES(IDIM) must be .GE. 4, but may be as +C large as the arrays COEF and WORK allow. +C The dimension is assumed to be NODES(NDIM). +C +C NOTE: The node grid is completely defined by +C the arguments XMIN, XMAX and NODES. +C The spacing of this grid in the IDIM +C coordinate direction is: +C +C DX(IDIM) = (XMAX(IDIM)-XMIN(IDIM)) / +C (NODES(IDIM)-1). +C +C A node in this grid may be indexed by +C an NDIM vector of integers +C (IN(1),...,IN(NDIM)) where +C 1 .LE. IN(IDIM) .LE. NODES(IDIM). +C The location of such a node may be +C represented by an NDIM vector +C (X(1),...,X(NDIM)) where +C X(IDIM) = XMIN(IDIM) + (IN(IDIM)-1) * +C DX(IDIM). +C +C XTRAP +C A parameter to control extrapolation to data +C sparse areas. The region described by XMIN +C and XMAX is divided into rectangles, the +C number of which is determined by NODES, and +C any rectangle containing a disproportionately +C small number of data points is considered to +C be data sparse (rectangle is used here to +C mean NDIM-dimensional rectangle). If XTRAP +C is nonzero the least squares problem is +C augmented with derivative constraints in the +C data sparse areas to prevent the matrix from +C becoming poorly conditioned. XTRAP serves as +C a weight for these constraints, and thus may +C be used to control smoothness in data sparse +C areas. Experience indicates that unity is a +C good first guess for this parameter. +C +C NOTE: If XTRAP is zero, substantial +C portions of the routine will be +C skipped, but a singular matrix +C can result if large portions of +C the region are without data. +C +C NCF +C The length of the array COEF in the calling +C program. If NCF is .LT. +C NODES(1)*...*NODES(NDIM), a fatal error is +C diagnosed. +C +C WORK +C A workspace array for solving the least +C squares matrix generated by this routine. +C Its required size is a function of the total +C number of nodes in the node grid. This +C total, NCOL = NODES(1)*...*NODES(NDIM), is +C also the number of columns in the least +C squares matrix. The length of the array WORK +C must equal or exceed NCOL*(NCOL+1). +C +C NWRK +C The length of the array WORK in the calling +C program. If +C NCOL = NODES(1)*...*NODES(NDIM) is the total +C number of nodes, then a fatal error is +C diagnosed if NWRK is less than +C NCOL*(NCOL+1). +C +C ON OUTPUT COEF +C The array of coefficients computed by this +C routine. Each coefficient corresponds to a +C particular basis function which in turn +C corresponds to a node in the node grid. This +C correspondence between the node grid and the +C array COEF is as if COEF were an +C NDIM-dimensional Fortran array with +C dimensions NODES(1),...,NODES(NDIM), i.e., to +C store the array linearly, the leftmost +C indices are incremented most frequently. +C Hence the length of the COEF array must equal +C or exceed the total number of nodes, which is +C NODES(1)*...*NODES(NDIM). The computed array +C COEF may be used with function SPLFE +C (or SPLDE) to evaluate the spline (or its +C derivatives) at an arbitrary point in NDIM +C*PL*ERROR* Comment line too long +C space. The dimension is assumed to be COEF(NCF). +C +C WORK +C The workspace containing intermediate +C calculations. It need not be saved. +C +C IERROR +C An error flag with the following meanings: +C 0 No error. +C 101 NDIM is .LT. 1 or is .GT. 4. +C 102 NODES(IDIM) is .LT. 4 fOR some IDIM. +C 103 XMIN(IDIM) = XMAX(IDIM) for some IDIM. +C 104 NCF (size of COEF) is +C .LT. NODES(1)*...*NODES(NDIM). +C 105 NDATA is .LT. 1. +C 106 NWRK (size of WORK) is too small. +C 107 SUPRLS failure (usually insufficient +C data) -- ordinarily occurs only if +C XTRAP is zero or WDATA contains all +C zeros. +C +C ALGORITHM An overdetermined system of linear equations +C is formed -- one equation for each data point +C plus equations for derivative constraints. +C This system is solved using subroutine SUPRLSD. +C +C ACCURACY If there is exactly one data point in the +C near vicinity of each node and no extra data, +C the resulting spline will agree with the +C data values to machine accuracy. However, if +C the problem is overdetermined or the sparse +C data option is utilized, the accuracy is hard +C to predict. Basically, smooth functions +C require fewer nodes than rough ones for the +C same accuracy. +C +C TIMING The execution time is roughly proportional +C to NDATA*NCOF**2 where NCOF = NODES(1)*...* +C NODES(NDIM). +C +C*********************************************************************** +C + SUBROUTINE SPLCCD(NDIM,XDATA,L1XDAT,YDATA,NDATA,XMIN,XMAX,NODES, + + XTRAP,COEF,NCF,WORK,NWRK,IERROR) + DOUBLE PRECISION XDATA + DOUBLE PRECISION YDATA + DOUBLE PRECISION XMIN + DOUBLE PRECISION XMAX + DOUBLE PRECISION XTRAP + DOUBLE PRECISION COEF + DOUBLE PRECISION WORK + DOUBLE PRECISION W + DIMENSION XDATA(L1XDAT,NDATA),YDATA(NDATA),XMIN(NDIM),XMAX(NDIM), + + NODES(NDIM),COEF(NCF),WORK(NWRK) + DIMENSION W(1) + SAVE +C + W(1) = -1.D0 + CALL SPLCWD(NDIM,XDATA,L1XDAT,YDATA,W,NDATA,XMIN,XMAX,NODES,XTRAP, + + COEF,NCF,WORK,NWRK,IERROR) +C + RETURN + END diff --git a/src/splcwd.f b/src/splcwd.f new file mode 100644 index 0000000..5b73079 --- /dev/null +++ b/src/splcwd.f @@ -0,0 +1,439 @@ +C +C Copyright (C) 2000 +C University Corporation for Atmospheric Research +C All Rights Reserved +C +C The use of this Software is governed by a License Agreement. +C + SUBROUTINE SPLCWD(NDIM,XDATA,L1XDAT,YDATA,WDATA,NDATA,XMIN,XMAX, + + NODES,XTRAP,COEF,NCF,WORK,NWRK,IERROR) + DOUBLE PRECISION XDATA + DOUBLE PRECISION YDATA + DOUBLE PRECISION WDATA + DOUBLE PRECISION XMIN + DOUBLE PRECISION XMAX + DOUBLE PRECISION XTRAP + DOUBLE PRECISION COEF + DOUBLE PRECISION WORK + DOUBLE PRECISION X + DOUBLE PRECISION DX + DOUBLE PRECISION DXIN + DOUBLE PRECISION SPCRIT + DOUBLE PRECISION XRNG + DOUBLE PRECISION SWGHT + DOUBLE PRECISION ROWWT + DOUBLE PRECISION RHS + DOUBLE PRECISION BASM + DOUBLE PRECISION RESERR + DOUBLE PRECISION TOTLWT + DOUBLE PRECISION BUMP + DOUBLE PRECISION WTPRRC + DOUBLE PRECISION EXPECT + DOUBLE PRECISION DCWGHT + DIMENSION XDATA(L1XDAT,NDATA),YDATA(NDATA),WDATA(NDATA), + + XMIN(NDIM),XMAX(NDIM),NODES(NDIM),COEF(NCF),WORK(NWRK) + DIMENSION X(4),NDERIV(4),IN(4),INMX(4) + COMMON /SPLCOMD/DX(4),DXIN(4),MDIM,IB(4),IBMN(4),IBMX(4) + SAVE +C +C The restriction that NDIM be less than are equat to 4 can be +C eliminated by increasing the above dimensions, but the required +C length of WORK becomes quite large. +C +C SPCRIT is used to determine data sparseness as follows - +C the weights assigned to all data points are totaled into the +C variable TOTLWT. (If no weights are entered, it is set to +C NDATA.) Each node of the node network is assigned a +C rectangle (in which it is contained) and the weights of all +C data points which fall in that rectangle are totaled. If that +C total is less than SPCRIT*EXPECT (EXPECT is defined below), +C then the node is ascertained to be in a data sparse location. +C EXPECT is that fraction of TOTLWT that would be expected by +C comparing the area of the rectangle with the total area under +C consideration. +C + DATA SPCRIT/.75D0/ +C + IERROR = 0 + MDIM = NDIM + IF (MDIM.LT.1 .OR. MDIM.GT.4) GO TO 127 + NCOL = 1 + DO 101 IDIM = 1,MDIM + NOD = NODES(IDIM) + IF (NOD.LT.4) GO TO 128 +C +C Number of columns in least squares matrix = number of coefficients = +C product of nodes over all dimensions. +C + NCOL = NCOL*NOD + XRNG = XMAX(IDIM) - XMIN(IDIM) + IF (XRNG.EQ.0.D0) GO TO 129 +C +C DX(IDIM) is the node spacing along the IDIM coordinate. +C + DX(IDIM) = XRNG/DBLE(NOD-1) + DXIN(IDIM) = 1.D0/DX(IDIM) + NDERIV(IDIM) = 0 + 101 CONTINUE + IF (NCOL.GT.NCF) GO TO 130 + NWRK1 = 1 + MDATA = NDATA + IF (MDATA.LT.1) GO TO 131 +C +C SWGHT is a local variable = XTRAP, and can be considered a smoothing +C weight for data sparse areas. If SWGHT .EQ. 0, no smoothing +C computations are performed. +C + SWGHT = XTRAP +C +C Set aside workspace for counting data points. +C + IF (SWGHT.NE.0.D0) NWRK1 = NCOL + 1 +C +C NWLFT is the length of the remaining workspace. +C + NWLFT = NWRK - NWRK1 + 1 + IF (NWLFT.LT.1) GO TO 132 + IROW = 0 +C +C ROWWT is used to weight rows of the least squares matrix. +C + ROWWT = 1.D0 +C +C Loop through all data points, computing a row for each. +C + DO 108 IDATA = 1,MDATA +C +C WDATA(1).LT.0 means weights have not been entered. In that case, +C ROWWT is left equal to 1. for all points. Otherwise ROWWT is +C equal to WDATA(IDATA). +C +C Every element of the row, as well as the corresponding right hand +C side, is multiplied by ROWWT. +C + IF (WDATA(1).LT.0.D0) GO TO 102 + ROWWT = WDATA(IDATA) +C +C Data points with 0 weight are ignored. +C + IF (ROWWT.EQ.0.D0) GO TO 108 + 102 IROW = IROW + 1 +C +C One row of the least squares matrix corresponds to each data +C point. The right hand for that row will correspond to the +C function value YDATA at that point. +C + RHS = ROWWT*YDATA(IDATA) + DO 103 IDIM = 1,MDIM + X(IDIM) = XDATA(IDIM,IDATA) + 103 CONTINUE +C +C The COEF array serves as a row of least squares matrix. +C Its value is zero except for columns corresponding to functions +C which are nonzero at X. +C + DO 104 ICOL = 1,NCOL + COEF(ICOL) = 0.D0 + 104 CONTINUE +C +C Compute the indices of basis functions which are nonzero at X. +C IBMN is in the range 0 to nodes-2 and IBMX is in range 1 +C to NODES-1. +C + DO 105 IDIM = 1,MDIM + NOD = NODES(IDIM) + IT = DXIN(IDIM)* (X(IDIM)-XMIN(IDIM)) + IBMN(IDIM) = MIN0(MAX0(IT-1,0),NOD-2) + IB(IDIM) = IBMN(IDIM) + IBMX(IDIM) = MAX0(MIN0(IT+2,NOD-1),1) + 105 CONTINUE +C +C Begining of basis index loop - traverse all indices corresponding +C to basis functions which are nonzero at X. The indices are in +C IB and are passed through common to BASCMP. +C + 106 CALL BASCMPD(X,NDERIV,XMIN,NODES,ICOL,BASM) +C +C BASCMP computes ICOL and BASM where BASM is the value at X of +C the N-dimensional basis function corresponding to column ICOL. +C + COEF(ICOL) = ROWWT*BASM +C +C Increment the basis indices. +C + DO 107 IDIM = 1,MDIM + IB(IDIM) = IB(IDIM) + 1 + IF (IB(IDIM).LE.IBMX(IDIM)) GO TO 106 + IB(IDIM) = IBMN(IDIM) + 107 CONTINUE +C +C End of basis index loop. +C +C +C Send a row of the least squares matrix to the reduction routine. +C + CALL SUPRLD(IROW,COEF,NCOL,RHS,WORK(NWRK1),NWLFT,COEF,RESERR, + + LSERR) + IF (LSERR.NE.0) GO TO 133 + 108 CONTINUE +C +C Row computations for all data points are now complete. +C +C If SWGHT.EQ.0, the least squares matrix is complete and no +C smoothing rows are computed. +C + IF (SWGHT.EQ.0.D0) GO TO 126 +C +C Initialize smoothing computations for data sparse areas. +C Derivative constraints will always have zero right hand side. +C + RHS = 0.D0 + NRECT = 1 +C +C Initialize the node indices and compute number of rectangles +C formed by the node network. +C + DO 109 IDIM = 1,MDIM + IN(IDIM) = 0 + INMX(IDIM) = NODES(IDIM) - 1 + NRECT = NRECT*INMX(IDIM) + 109 CONTINUE +C +C Every node is assigned an element of the workspace (set aside +C previously) in which data points are counted. +C + DO 110 IIN = 1,NCOL + WORK(IIN) = 0.D0 + 110 CONTINUE +C +C Assign each data point to a node, total the assignments for +C each node, and save in the workspace. +C + TOTLWT = 0.D0 + DO 112 IDATA = 1,MDATA +C +C BUMP is the weight associated with the data point. +C + BUMP = 1.D0 + IF (WDATA(1).GE.0.D0) BUMP = WDATA(IDATA) + IF (BUMP.EQ.0.D0) GO TO 112 +C +C Find the nearest node. +C + IIN = 0 + DO 111 IDIMC = 1,MDIM + IDIM = MDIM + 1 - IDIMC + INIDIM = INT(DXIN(IDIM)* (XDATA(IDIM,IDATA)-XMIN(IDIM))+ + + .5D0) +C +C Points not in range (+ or - 1/2 node spacing) are not counted. +C + IF (INIDIM.LT.0 .OR. INIDIM.GT.INMX(IDIM)) GO TO 112 +C +C Compute linear address of node in workspace by Horner's method. +C + IIN = (INMX(IDIM)+1)*IIN + INIDIM + 111 CONTINUE +C +C Bump counter for that node. +C + WORK(IIN+1) = WORK(IIN+1) + BUMP + TOTLWT = TOTLWT + BUMP + 112 CONTINUE +C +C Compute the expected weight per rectangle. +C + WTPRRC = TOTLWT/DBLE(NRECT) +C +C IN contains indices of the node (previously initialized). +C IIN will be the linear address of the node in the workspace. +C + IIN = 0 +C +C Loop through all nodes, computing derivative constraint rows +C for those in data sparse locations. +C +C Begining of node index loop - traverse all node indices. +C The indices are in IN. +C + 113 IIN = IIN + 1 + EXPECT = WTPRRC +C +C Rectangles at edge of network are smaller and hence less weight +C should be expected. +C + DO 114 IDIM = 1,MDIM + IF (IN(IDIM).EQ.0 .OR. IN(IDIM).EQ.INMX(IDIM)) EXPECT = .5D0* + + EXPECT + 114 CONTINUE +C +C The expected weight minus the actual weight serves to define +C data sparseness and is also used to weight the derivative +C constraint rows. +C +C There is no constraint if not data sparse. +C + IF (WORK(IIN).GE.SPCRIT*EXPECT) GO TO 124 + DCWGHT = EXPECT - WORK(IIN) + DO 115 IDIM = 1,MDIM + INIDIM = IN(IDIM) +C +C Compute the location of the node. +C + X(IDIM) = XMIN(IDIM) + DBLE(INIDIM)*DX(IDIM) +C +C Compute the indices of the basis functions which are non-zero +C at the node. +C + IBMN(IDIM) = INIDIM - 1 + IBMX(IDIM) = INIDIM + 1 +C +C Distinguish the boundaries. +C + IF (INIDIM.EQ.0) IBMN(IDIM) = 0 + IF (INIDIM.EQ.INMX(IDIM)) IBMX(IDIM) = INMX(IDIM) +C +C Initialize the basis indices. +C + IB(IDIM) = IBMN(IDIM) + 115 CONTINUE +C +C Multiply by the extrapolation parameter (this acts as a +C smoothing weight). +C + DCWGHT = SWGHT*DCWGHT +C +C The COEF array serves as a row of the least squares matrix. +C Its value is zero except for columns corresponding to functions +C which are non-zero at the node. +C + DO 116 ICOL = 1,NCOL + COEF(ICOL) = 0.D0 + 116 CONTINUE +C +C The 2nd derivative of a function of MDIM variables may be thought +C of as a symmetric MDIM x MDIM matrix of 2nd order partial +C derivatives. Traverse the upper triangle of this matrix and, +C for each element, compute a row of the least squares matrix. +C + DO 123 IDM = 1,MDIM + DO 122 JDM = IDM,MDIM + DO 117 IDIM = 1,MDIM + NDERIV(IDIM) = 0 + 117 CONTINUE +C +C Off-diagonal elements appear twice by symmetry, so the corresponding +C row is weighted by a factor of 2. +C + ROWWT = 2.D0*DCWGHT + IF (JDM.NE.IDM) GO TO 118 +C +C Diagonal. +C + ROWWT = DCWGHT + NDERIV(JDM) = 2 + IF (IN(IDM).NE.0 .AND. IN(IDM).NE.INMX(IDM)) GO TO 119 +C +C Node is at boundary. +C +C Normal 2nd derivative constraint at boundary is not appropriate for +C natural splines (2nd derivative 0 by definition). Substitute +C a 1st derivative constraint. +C + 118 NDERIV(IDM) = 1 + NDERIV(JDM) = 1 + 119 IROW = IROW + 1 +C +C Begining of basis index loop - traverse all indices corresponding +C to basis functions which are non-zero at X. +C The indices are in IB and are passed through common to BASCMP. +C + 120 CALL BASCMPD(X,NDERIV,XMIN,NODES,ICOL,BASM) +C +C BASCMP computes ICOL and BASM where BASM is the value at X of the +C N-dimensional basis function corresponding to column ICOL. +C + COEF(ICOL) = ROWWT*BASM +C +C Increment the basis indices. +C + DO 121 IDIM = 1,MDIM + IB(IDIM) = IB(IDIM) + 1 + IF (IB(IDIM).LE.IBMX(IDIM)) GO TO 120 + IB(IDIM) = IBMN(IDIM) + 121 CONTINUE +C +C End of basis index loop. +C +C Send row of least squares matrix to reduction routine. +C + CALL SUPRLD(IROW,COEF,NCOL,RHS,WORK(NWRK1),NWLFT,COEF, + + RESERR,LSERR) + IF (LSERR.NE.0) GO TO 133 + 122 CONTINUE + 123 CONTINUE +C +C Increment node indices. +C + 124 DO 125 IDIM = 1,MDIM + IN(IDIM) = IN(IDIM) + 1 + IF (IN(IDIM).LE.INMX(IDIM)) GO TO 113 + IN(IDIM) = 0 + 125 CONTINUE +C +C End of node index loop. +C +C Call for least squares solution in COEF array. +C + 126 IROW = 0 + CALL SUPRLD(IROW,COEF,NCOL,RHS,WORK(NWRK1),NWLFT,COEF,RESERR, + + LSERR) + IF (LSERR.NE.0) GO TO 133 + RETURN +C +C Error section +C + 127 CONTINUE + IERROR = 101 + CALL CFAERR(IERROR, + + ' SPLCCD or SPLCWD - NDIM is less than 1 or is greater than 4' + + ,60) + GO TO 134 + 128 CONTINUE + IERROR = 102 + CALL CFAERR(IERROR, + + ' SPLCCD or SPLCWD - NODES(IDIM) is less than 4 for some IDIM' + + ,60) + GO TO 134 + 129 CONTINUE + IERROR = 103 + CALL CFAERR(IERROR, + + ' SPLCCD or SPLCWD - XMIN(IDIM) equals XMAX(IDIM) for some IDIM' + + ,60) + GO TO 134 + 130 CONTINUE + IERROR = 104 + CALL CFAERR(IERROR, + + ' SPLCCD or SPLCWD - NCF (size of COEF) is too small ' + + ,60) + GO TO 134 + 131 CONTINUE + IERROR = 105 + CALL CFAERR(IERROR, + + ' SPLCCD or SPLCWD - Ndata Is less than 1 ' + + ,60) + GO TO 134 + 132 CONTINUE + IERROR = 106 + CALL CFAERR(IERROR, + + ' SPLCCD or SPLCWD - NWRK (size of WORK) is too small ' + + ,60) + GO TO 134 + 133 CONTINUE + IERROR = 107 + CALL CFAERR(IERROR, + +' SPLCCD or SPLCWD - SUPRLS failure (this usually indicates insuff + +icient input data',80) +C + 134 RETURN + END diff --git a/src/splded.f b/src/splded.f new file mode 100644 index 0000000..678b296 --- /dev/null +++ b/src/splded.f @@ -0,0 +1,116 @@ +C +C Copyright (C) 2000 +C University Corporation for Atmospheric Research +C All Rights Reserved +C +C The use of this Software is governed by a License Agreement. +C + FUNCTION SPLDED(NDIM,X,NDERIV,COEF,XMIN,XMAX,NODES,IERROR) + DOUBLE PRECISION SPLDED + DOUBLE PRECISION X + DOUBLE PRECISION COEF + DOUBLE PRECISION XMIN + DOUBLE PRECISION XMAX + DOUBLE PRECISION DX + DOUBLE PRECISION DXIN + DOUBLE PRECISION XRNG + DOUBLE PRECISION SUM + DOUBLE PRECISION BASM + DIMENSION X(NDIM),NDERIV(NDIM),COEF(*),XMIN(NDIM),XMAX(NDIM), + + NODES(NDIM) + COMMON /SPLCOMD/DX(4),DXIN(4),MDIM,IB(4),IBMN(4),IBMX(4) + SAVE +C +C The restriction for NDIM to be .LE. 4 can be eliminated by increasing +C the above dimensions. +C + IERROR = 0 + MDIM = NDIM + IF (MDIM.LT.1 .OR. MDIM.GT.4) GO TO 105 + IIBMX = 1 + DO 101 IDIM = 1,MDIM + NOD = NODES(IDIM) + IF (NOD.LT.4) GO TO 106 + XRNG = XMAX(IDIM) - XMIN(IDIM) + IF (XRNG.EQ.0.D0) GO TO 107 + IF (NDERIV(IDIM).LT.0 .OR. NDERIV(IDIM).GT.2) GO TO 108 +C +C DX(IDIM) is the node spacing along the IDIM coordinate. +C + DX(IDIM) = XRNG/DBLE(NOD-1) + DXIN(IDIM) = 1.D0/DX(IDIM) +C +C Compute indices of basis functions which are nonzero at X. +C + IT = DXIN(IDIM)* (X(IDIM)-XMIN(IDIM)) +C +C IBMN must be in the range 0 to NODES-2. +C + IBMN(IDIM) = MIN0(MAX0(IT-1,0),NOD-2) +C +C IBMX must be in the range 1 to NODES-1. +C + IBMX(IDIM) = MAX0(MIN0(IT+2,NOD-1),1) + IIBMX = IIBMX* (IBMX(IDIM)-IBMN(IDIM)+1) + IB(IDIM) = IBMN(IDIM) + 101 CONTINUE +C + SUM = 0.D0 + IIB = 0 +C +C Begining of basis index loop - traverse all indices corresponding +C to basis functions which are nonzero at X. +C + 102 IIB = IIB + 1 +C +C The indices are in IB and are passed through common to BASCMP. +C + CALL BASCMPD(X,NDERIV,XMIN,NODES,ICOF,BASM) +C +C BASCMP computes ICOF and BASM where BASM is the value at X of the +C N-dimensional basis function corresponding to COEF(ICOF). +C + SUM = SUM + COEF(ICOF)*BASM + IF (IIB.GE.IIBMX) GO TO 104 +C +C Increment the basis indices. +C + DO 103 IDIM = 1,MDIM + IB(IDIM) = IB(IDIM) + 1 + IF (IB(IDIM).LE.IBMX(IDIM)) GO TO 102 + IB(IDIM) = IBMN(IDIM) + 103 CONTINUE +C +C End of basis index loop. +C + 104 SPLDED = SUM + RETURN +C +C Errors. +C + 105 CONTINUE + IERROR = 101 + CALL CFAERR(IERROR, + + ' SPLFED or SPLDED - NDIM is less than 1 or greater than 4 ' + + ,60) + GO TO 109 + 106 CONTINUE + IERROR = 102 + CALL CFAERR(IERROR, + + ' SPLFED or SPLDED - NODES(IDIM) is less than 4for some IDIM' + + ,60) + GO TO 109 + 107 CONTINUE + IERROR = 103 + CALL CFAERR(IERROR, + + ' SPLFED or SPLDED - XMIN(IDIM) = XMAX(IDIM) for some IDIM ' + + ,60) + GO TO 109 + 108 CONTINUE + IERROR = 104 + CALL CFAERR(IERROR, + +' SPLDED - NDERIV(IDIM) IS less than 0 or greater than 2 for some + +IDIM ',70) +C + 109 STOP + END diff --git a/src/splfed.f b/src/splfed.f new file mode 100644 index 0000000..26250a9 --- /dev/null +++ b/src/splfed.f @@ -0,0 +1,27 @@ +C +C Copyright (C) 2000 +C University Corporation for Atmospheric Research +C All Rights Reserved +C +C The use of this Software is governed by a License Agreement. +C + FUNCTION SPLFED(NDIM,X,COEF,XMIN,XMAX,NODES,IERROR) + DOUBLE PRECISION SPLFED + DOUBLE PRECISION X + DOUBLE PRECISION COEF + DOUBLE PRECISION XMIN + DOUBLE PRECISION XMAX + DOUBLE PRECISION SPLDED + DIMENSION X(NDIM),COEF(*),XMIN(NDIM),XMAX(NDIM),NODES(NDIM) + DIMENSION NDERIV(4) + SAVE +C + DATA NDERIV(1),NDERIV(2),NDERIV(3),NDERIV(4)/0,0,0,0/ +C +C The restriction for NDIM to be .LE. 4 can be eliminated by +C increasing the above dimension and those in SPLDED. +C + SPLFED = SPLDED(NDIM,X,NDERIV,COEF,XMIN,XMAX,NODES,IERROR) +C + RETURN + END diff --git a/src/suprld.f b/src/suprld.f new file mode 100644 index 0000000..7ac2f79 --- /dev/null +++ b/src/suprld.f @@ -0,0 +1,288 @@ +C +C Copyright (C) 2000 +C University Corporation for Atmospheric Research +C All Rights Reserved +C +C The use of this Software is governed by a License Agreement. +C + SUBROUTINE SUPRLD(I,ROWI,N,BI,A,NN,SOLN,ERR,IER) + DOUBLE PRECISION ROWI + DOUBLE PRECISION BI + DOUBLE PRECISION A + DOUBLE PRECISION SOLN + DOUBLE PRECISION ERR + DOUBLE PRECISION ERRSUM + DOUBLE PRECISION S + DOUBLE PRECISION TEMP + DOUBLE PRECISION TEMP1 + DOUBLE PRECISION CN + DOUBLE PRECISION SN + DIMENSION ROWI(N),A(NN),SOLN(N) + SAVE +C + IER = 0 + IF (I.GT.1) GO TO 101 +C +C Routine entered with I.LE.0 means complete the reduction and store +C the solution in SOLN. +C + IF (I.LE.0) GO TO 125 +C +C Set up quantities on first call. +C + IOLD = 0 + NP1 = N + 1 +C +C Compute how many rows can be input now. +C + L = NN/NP1 + ILAST = 0 + IL1 = 0 + K = 0 + K1 = 0 + ERRSUM = 0.D0 + NREQ = ((N+5)*N+2)/2 +C +C Error exit if insufficient scratch storage provided. +C + IF (NN.GE.NREQ) GO TO 101 + IER = 32 + CALL CFAERR(IER, + +' SUPRLD - insufficient scratch storage provided. at least ((N+5)* + +N+2)/2 locations needed',88) + RETURN +C +C Store the row in the scratch storage. +C + 101 CONTINUE +C +C Error exit if (I-IOLD).NE.1. +C + IF ((I-IOLD).EQ.1) GO TO 102 + IER = 35 + CALL CFAERR(IER,' SUPRLD - values of I not in sequence',37) + RETURN +C + 102 CONTINUE + IOLD = I + DO 103 J = 1,N + ILJ = ILAST + J + A(ILJ) = ROWI(J) + 103 CONTINUE + ILNP = ILAST + NP1 + A(ILNP) = BI + ILAST = ILAST + NP1 + ISAV = I + IF (I.LT.L) RETURN + 104 CONTINUE + IF (K.EQ.0) GO TO 115 + K1 = MIN0(K,N) + IDIAG = -NP1 + IF (L-K.EQ.1) GO TO 110 +C +C Apply householder transformations to zero out new rows. +C + DO 109 J = 1,K1 + IDIAG = IDIAG + (NP1-J+2) + I1 = IL1 + J + I2 = I1 + NP1* (L-K-1) + S = A(IDIAG)*A(IDIAG) + DO 105 II = I1,I2,NP1 + S = S + A(II)*A(II) + 105 CONTINUE + IF (S.EQ.0.D0) GO TO 109 + TEMP = A(IDIAG) + A(IDIAG) = SQRT(S) + IF (TEMP.GT.0.D0) A(IDIAG) = -A(IDIAG) + TEMP = TEMP - A(IDIAG) + TEMP1 = 1.D0/ (TEMP*A(IDIAG)) + JP1 = J + 1 + DO 108 J1 = JP1,NP1 + JDEL = J1 - J + IDJ = IDIAG + JDEL + S = TEMP*A(IDJ) + DO 106 II = I1,I2,NP1 + IIJD = II + JDEL + S = S + A(II)*A(IIJD) + 106 CONTINUE + S = S*TEMP1 + A(IDJ) = A(IDJ) + S*TEMP + DO 107 II = I1,I2,NP1 + IIJD = II + JDEL + A(IIJD) = A(IIJD) + S*A(II) + 107 CONTINUE + 108 CONTINUE + 109 CONTINUE + GO TO 113 +C +C Apply rotations to zero out the single new row. +C + 110 DO 112 J = 1,K1 + IDIAG = IDIAG + (NP1-J+2) + I1 = IL1 + J + IF (ABS(A(I1)).LE.1.D-18) THEN + S = SQRT(A(IDIAG)*A(IDIAG)) + ELSE IF (ABS(A(IDIAG)).LT.1.D-18) THEN + S = SQRT(A(I1)*A(I1)) + ELSE + S = SQRT(A(IDIAG)*A(IDIAG)+A(I1)*A(I1)) + END IF + IF (S.EQ.0.D0) GO TO 112 + TEMP = A(IDIAG) + A(IDIAG) = S + S = 1.D0/S + CN = TEMP*S + SN = A(I1)*S + JP1 = J + 1 + DO 111 J1 = JP1,NP1 + JDEL = J1 - J + IDJ = IDIAG + JDEL + TEMP = A(IDJ) + I1JD = I1 + JDEL + A(IDJ) = CN*TEMP + SN*A(I1JD) + A(I1JD) = -SN*TEMP + CN*A(I1JD) + 111 CONTINUE + 112 CONTINUE + 113 IF (K.LT.N) GO TO 115 + LMKM1 = L - K +C +C Accumulate residual sum of squares. +C + DO 114 II = 1,LMKM1 + ILNP = IL1 + II*NP1 + ERRSUM = ERRSUM + A(ILNP)*A(ILNP) + 114 CONTINUE + IF (I.LE.0) GO TO 127 + K = L + ILAST = IL1 +C +C Determine how many new rows may be input on next iteration. +C + L = K + (NN-ILAST)/NP1 + RETURN + 115 K11 = K1 + 1 + K1 = MIN0(L,N) + IF (L-K.EQ.1) GO TO 122 + K1M1 = K1 - 1 + IF (L.GT.N) K1M1 = N + I1 = IL1 + K11 - NP1 - 1 +C +C Perform householder transformations to reduce rows to upper +C triangular form. +C + DO 120 J = K11,K1M1 + I1 = I1 + (NP1+1) + I2 = I1 + (L-J)*NP1 + S = 0.D0 + DO 116 II = I1,I2,NP1 + S = S + A(II)*A(II) + 116 CONTINUE + IF (S.EQ.0.D0) GO TO 120 + TEMP = A(I1) + A(I1) = SQRT(S) + IF (TEMP.GT.0.D0) A(I1) = -A(I1) + TEMP = TEMP - A(I1) + TEMP1 = 1.D0/ (TEMP*A(I1)) + JP1 = J + 1 + I11 = I1 + NP1 + DO 119 J1 = JP1,NP1 + JDEL = J1 - J + I1JD = I1 + JDEL + S = TEMP*A(I1JD) + DO 117 II = I11,I2,NP1 + IIJD = II + JDEL + S = S + A(II)*A(IIJD) + 117 CONTINUE + S = S*TEMP1 + I1JD = I1 + JDEL + A(I1JD) = A(I1JD) + S*TEMP + DO 118 II = I11,I2,NP1 + IIJD = II + JDEL + A(IIJD) = A(IIJD) + S*A(II) + 118 CONTINUE + 119 CONTINUE + 120 CONTINUE + IF (L.LE.N) GO TO 122 + NP1MK = NP1 - K + LMK = L - K +C +C Accumulate residual sum of squares. +C + DO 121 II = NP1MK,LMK + ILNP = IL1 + II*NP1 + ERRSUM = ERRSUM + A(ILNP)*A(ILNP) + 121 CONTINUE + 122 IMOV = 0 + I1 = IL1 + K11 - NP1 - 1 +C +C Squeeze the unnecessary elements out of scratch storage to +C allow space for more rows. +C + DO 124 II = K11,K1 + IMOV = IMOV + (II-1) + I1 = I1 + NP1 + 1 + I2 = I1 + NP1 - II + DO 123 III = I1,I2 + IIIM = III - IMOV + A(IIIM) = A(III) + 123 CONTINUE + 124 CONTINUE + ILAST = I2 - IMOV + IL1 = ILAST + IF (I.LE.0) GO TO 127 + K = L +C +C Determine how many new rows may be input on next iteration. +C + L = K + (NN-ILAST)/NP1 + RETURN +C +C Complete reduction and store solution in SOLN. +C + 125 L = ISAV +C +C Error exit if L less than N. +C + IF (L.GE.N) GO TO 126 + IER = 33 + CALL CFAERR(IER,' SUPRLD - array has too few rows.',33) + RETURN + 126 CONTINUE +C +C K.NE.ISAV means further reduction needed. +C + IF (K.NE.ISAV) GO TO 104 + 127 ILAST = (NP1* (NP1+1))/2 - 1 + IF (A(ILAST-1).EQ.0.D0) GO TO 130 +C +C Solve triangular system into ROWI. +C + SOLN(N) = A(ILAST)/A(ILAST-1) + DO 129 II = 2,N + IIM1 = II - 1 + ILAST = ILAST - II + S = A(ILAST) + DO 128 K = 1,IIM1 + ILK = ILAST - K + NPK = NP1 - K + S = S - A(ILK)*SOLN(NPK) + 128 CONTINUE + ILII = ILAST - II + IF (A(ILII).EQ.0.D0) GO TO 130 + NPII = NP1 - II + SOLN(NPII) = S/A(ILII) + 129 CONTINUE +C +C Store residual norm. +C + ERR = SQRT(ERRSUM) + RETURN +C +C Error return if system is singular. +C + 130 CONTINUE + IER = 34 + CALL CFAERR(IER,' SUPRLD - system is singular.',29) + RETURN +C + END