[newfield_revision Patch] Add Linear Algebra to Statigraphic Flow
Jeffrey Oldham
oldham at codesourcery.com
Thu Aug 16 20:31:15 UTC 2001
This patch to the (Chevron^r) statigraphic flow program adds the
linear algebra computation. We use the open-source beta Template
Numerical Toolkit package from the NIST. The linear algebra
computation needs to be moved to a stencil. We are currently working
on that.
2001-08-16 Jeffrey D. Oldham <oldham at codesourcery.com>
* StatigraphicFlow.cpp: Add header files for linear algebra
computations. Remove answered question and finished functions.
Add sum() to unfinished work list. Update explanations with
revised replicate() information and new information for
nearestNeighbors(Centering, FieldOffset, Centering) and sum().
(findIndex): New function.
(main): Fix creation of spoke centering. Add creation of disspoke
centering. Add directedPermeability and faceDistance fields.
Replace linear algebra pseudocode with explicit computation.
Revise computation of totalFlux.
* makefile: (RULE_INCLUDES): Add support for linear alegebra
library.
* tnt/A10x10.dat: New file from linear algebra package.
* tnt/AB10x10.dat: Likewise.
* tnt/Makefile: Likewise.
* tnt/SPD10x10.dat: Likewise.
* tnt/cholesky.h: Likewise.
* tnt/cmat.h: Likewise.
* tnt/fchol.cc: Likewise.
* tnt/fcscmat.h: Likewise.
* tnt/fmat.h: Likewise.
* tnt/fortran.h: Likewise.
* tnt/fspvec.h: Likewise.
* tnt/index.h: Likewise.
* tnt/lapack.h: Likewise.
* tnt/lu-C.cc: Likewise.
* tnt/lu.cc: Likewise.
* tnt/lu.h: Likewise.
* tnt/matrix.dat: Likewise.
* tnt/qr.cc: Likewise.
* tnt/qr.h: Likewise.
* tnt/region1d.h: Likewise.
* tnt/region2d.h: Likewise.
* tnt/stopwatch.h: Likewise.
* tnt/subscript.h: Likewise.
* tnt/tnt.h: Likewise.
* tnt/tnt.h.patch: Likewise.
* tnt/tntmath.h: Likewise.
* tnt/tntreqs.h: Likewise.
* tnt/transv.h: Likewise.
* tnt/triang.h: Likewise.
* tnt/trisolve.h: Likewise.
* tnt/vec.h: Likewise.
* tnt/vec.h.patch: Likewise.
* tnt/vecadaptor.h: Likewise.
* tnt/version.h: Likewise.
Applied to newfield_revision branch of examples/NewField/StatigraphicFlow/
Approved by Stephen Smith
Not tested. It does not currently cleanly compile.
Thanks,
Jeffrey D. Oldham
oldham at codesourcery.com
-------------- next part --------------
Index: StatigraphicFlow.cpp
===================================================================
RCS file: /home/pooma/Repository/r2/examples/NewField/StatigraphicFlow/Attic/StatigraphicFlow.cpp,v
retrieving revision 1.1.2.1
diff -c -p -r1.1.2.1 StatigraphicFlow.cpp
*** StatigraphicFlow.cpp 2001/08/04 02:19:06 1.1.2.1
--- StatigraphicFlow.cpp 2001/08/16 19:02:55
***************
*** 7,12 ****
--- 7,17 ----
#include <iostream>
#include <cstdlib>
#include "Pooma/NewFields.h"
+ // linear algebra files:
+ #include "tnt/tnt.h"
+ #include "tnt/vec.h"
+ #include "tnt/cmat.h"
+ #include "tnt/lu.h"
// This program implements "Implementation of a Flux-Continuous Fnite
// Difference Method for Stratigraphic, Hexahedron Grids," by
***************
*** 28,49 ****
// o. I omitted a separate coordinates field, presumably updated each
// iteration, in favor of using the mesh. Since I do not know how
// the coordinates are updated, I omitted updating the mesh.
- // o. Is it important to flesh out the linear algebra solution? We
- // might learn something about field syntax, but it will also take
- // time for me to determine the correct operands.
// o. Creating non-canonical edge and face centerings requires
// dimension-dependent code. Is this acceptable?
/** UNFINISHED WORK **/
- // o replicate(field, std::vector<FieldOffsetList>)
// o meshLayout.unitCoordinateNormals()
// o field.mesh()
// o field.mesh().normals()
// o field.mesh().normals().signedMagnitude()
-
/** EXPLANATIONS **/
// o Centering<Dim> canonicalCentering<Dim>(CellType, Continuous):
--- 33,50 ----
// o. I omitted a separate coordinates field, presumably updated each
// iteration, in favor of using the mesh. Since I do not know how
// the coordinates are updated, I omitted updating the mesh.
// o. Creating non-canonical edge and face centerings requires
// dimension-dependent code. Is this acceptable?
/** UNFINISHED WORK **/
// o meshLayout.unitCoordinateNormals()
// o field.mesh()
// o field.mesh().normals()
// o field.mesh().normals().signedMagnitude()
+ // o sum(field, vector<FieldOffsetList>, centering)
/** EXPLANATIONS **/
// o Centering<Dim> canonicalCentering<Dim>(CellType, Continuous):
***************
*** 62,72 ****
// o The Chevron algorithm first solves a linear program. I have
// omitted since computation since it does not illustrate field
// computations.
! // o replicate(field, std::vector<FieldOffsetList>): This function,
! // syntactic sugar for a nearest neighbors computation, copies the
! // field values to the positions indicated by the
! // std::vector<FieldOffsetList>. Each field value is copied to one
! // or more values. replicate() could be replaced by sum(), but the
// latter function has an unnecessary loop since each output value
// equals one input value.
// o nearestNeighbors(inputCentering, outputCentering): This function
--- 63,73 ----
// o The Chevron algorithm first solves a linear program. I have
// omitted since computation since it does not illustrate field
// computations.
! // o replicate(field, std::vector<FieldOffsetList>, centering): This
! // function, syntactic sugar for a nearest neighbors computation,
! // copies the field values to the positions indicated by the
! // std::vector<FieldOffsetList>. Each field value is copied to one or
! // more values. replicate() could be replaced by sum(), but the
// latter function has an unnecessary loop since each output value
// equals one input value.
// o nearestNeighbors(inputCentering, outputCentering): This function
***************
*** 75,80 ****
--- 76,87 ----
// value, the closest input values, wrt Manhattan distance, are
// returned. Eventually, these may be pre-computed or cached to
// reduce running time.
+ // o nearestNeighbors(input centering, FieldOffset, offset's centering
+ // [, bool]): This function returns a FieldOffsetList with entries for
+ // the closest input values, wrt Manhattan distance, to the value
+ // specified by the FieldOffset and the offset centering. The
+ // optional fourth parameter indicates only values from the FieldOffset's
+ // cell should be returned.
// o meshLayout.unitCoordinateNormals(): This returns a discontinuous
// face-centered field with unit-length normals all pointing in
// positive directions.
***************
*** 87,94 ****
// equalling the face's area/volume and sign equalling whether the
// face's normal is in a positive direction, e.g., the positive
// x-direction vs. the negative x-direction.
! // o sum(field, FieldOffsetList): this parallel-data statement adds
! // the values indicated in the FieldOffsetList to form each output value
/** DESIGN DECISIONS **/
--- 94,102 ----
// equalling the face's area/volume and sign equalling whether the
// face's normal is in a positive direction, e.g., the positive
// x-direction vs. the negative x-direction.
! // o sum(field, vector<FieldOffsetList>, centering): this
! // parallel-data statement adds the values indicated in the
! // FieldOffsetList to form each output value
/** DESIGN DECISIONS **/
***************
*** 98,103 ****
--- 106,129 ----
/** THE PROGRAM **/
+
+ // Return the index of the specified field offset in the given list.
+ // Return a negative number if not found.
+
+ template <int Dim>
+ inline int
+ findIndex(const FieldOffsetList<Dim> &vec,
+ const FieldOffset<Dim> &fo)
+ {
+ int indx;
+ for (indx = vec.size()-1;
+ indx >= 0 && vec[indx] != fo;
+ --indx)
+ ;
+ return indx;
+ }
+
+
int main(int argc, char *argv[])
{
// Set up the Pooma library.
*************** int main(int argc, char *argv[])
*** 158,172 ****
position(zeroFace) = 0.0;
position(1-zeroFace) = 0.25; spoke.addValue(orientation, position);
position(1-zeroFace) = 0.75; spoke.addValue(orientation, position);
- position(zeroFace) = 1.0;
- position(1-zeroFace) = 0.25; spoke.addValue(orientation, position);
- position(1-zeroFace) = 0.75; spoke.addValue(orientation, position);
}
Fields_t spokeFlux (spoke, meshLayout, origin, spacings);
// Face-centered.
Centering<Dim> disFace = canonicalCentering<Dim>(FaceType, Discontinuous);
/* INITIALIZATION */
--- 184,211 ----
position(zeroFace) = 0.0;
position(1-zeroFace) = 0.25; spoke.addValue(orientation, position);
position(1-zeroFace) = 0.75; spoke.addValue(orientation, position);
}
Fields_t spokeFlux (spoke, meshLayout, origin, spacings);
+ Centering<Dim> disspoke(FaceType, Discontinuous);
+ // NOTE: This code is not dimension-independent.
+ for (int zeroFace = 0; zeroFace < 2; ++zeroFace) {
+ orientation = 1; orientation[zeroFace] = 0;
+ position(zeroFace) = 0.0;
+ position(1-zeroFace) = 0.25; disspoke.addValue(orientation, position);
+ position(1-zeroFace) = 0.75; disspoke.addValue(orientation, position);
+ position(zeroFace) = 1.0;
+ position(1-zeroFace) = 0.25; disspoke.addValue(orientation, position);
+ position(1-zeroFace) = 0.75; disspoke.addValue(orientation, position);
+ }
+
// Face-centered.
Centering<Dim> disFace = canonicalCentering<Dim>(FaceType, Discontinuous);
+ Fieldv_t directedPermeability (disFace, meshLayout, origin, spacings);
+ // \gamma_{i,j} = K_i^t \dot \hat{n}_j
+ Fields_t faceDistance (disFace, meshLayout, origin, spacings);
+ // distance from cell center to face center
/* INITIALIZATION */
*************** int main(int argc, char *argv[])
*** 180,210 ****
/* COMPUTATION */
- #ifdef PSEUDOCODE
// Compute pressureGradients by simultaneously solving several
// linear equations. The operands have different centerings.
! // FIXME
! pressureGradients =
! linearAlgebra<2>(pressure /* cell-centered */,
! /* Interpolate from vertex-centered to cell-centered: */
! interpolate<Cell,Vertex>(coordinates),
! permeability /* cell-centered */,
! normals /* face-centered */);
#endif // PSEUDOCODE
! // Compute the spoke fluxes.
! // We must multiply three quantities, each with a different
! // centering, to yield values at a fourth-centering. permeability
! // is cell-centered. pressureGradient is subcell-centered. The
! // normals are face-centered. The product is spoke-centered.
spokeFlux =
! dot(replicate(dot(replicate(permeability, nearestNeighbors(cell, subcell)),
! pressureGradient),
! nearestNeighbors(subcell, spoke)),
! replicate(meshLayout.unitCoordinateNormals(),
! nearestNeighbors(disFace, spoke)));
// Sum the spoke fluxes into a cell flux.
--- 219,339 ----
/* COMPUTATION */
// Compute pressureGradients by simultaneously solving several
// linear equations. The operands have different centerings.
!
! // \Gamma is used in the flux continuity equations.
! directedPermeability =
! dot(replicate(permeability,
! nearestNeighbors(permeability.centering(),
! directedPermeability.centering()),
! directedPermeability.centering()),
! meshLayout.unitCoordinateNormals());
!
! #ifdef PSEUDOCODE
! // These distances are used in the pressure continuity equations.
! faceDistance = face-centered-positions -
! replicate(cell-centered-positions,
! nearestNeighbors(cell, faceDistance.centering()),
! faceDistance.centering());
#endif // PSEUDOCODE
+
+ const int nuRows = (1 << Dim) * Dim;
+ TNT::Matrix<double> A(nuRows, nuRows, 0.0);
+ TNT::Vector<double> rhs(nuRows, 0.0);
+ TNT::Vector<TNT::Subscript> ipiv; // ignored
+
+ // FIXME: Move this code to a stencil so it can be applied across
+ // the entire grid.
+
+ // Assign values to the matrix A and vector rhs.
+
+ const Centering<Dim> vert = canonicalCentering<Dim>(VertexType, Continuous);
+ PInsist(vert.size() == 1, "The vertex centering has too many values.");
+
+ const FieldOffsetList<Dim> gradients =
+ nearestNeighbors(pressureGradient.centering(),
+ FieldOffset<Dim>(Loc<Dim>(0)) /* cell origin */, vert);
+ // gradients's order of pressure gradients will be used for the
+ // matrix and rhs.
+
+ const FieldOffsetList<Dim> fluxPoints =
+ nearestNeighbors(spokeFlux.centering(),
+ FieldOffset<Dim>(Loc<Dim>(0)) /* cell origin */, vert);
+ // fluxPoints has locations for all faces incident to the vertex.
+ const int size = fluxPoints.size();
+ PAssert(gradients.size() * Dim == size * 2);
+
+ const std::vector<FieldOffsetList<Dim> > disFluxPoints =
+ nearestNeighbors(disspoke, fluxPoints, spokeFlux.centering());
+ // For every i, disFluxPoints[i] is two discontinuous positions on
+ // "either side" of the face represented by fluxPoints[i].
+
+ for (int faceIndex = size-1; faceIndex >= 0; --faceIndex) {
+ // Work on the "positive" side of the face.
+ FieldOffset<Dim> fo = disFluxPoints[faceIndex][0];
+ int columnNu =
+ findIndex(gradients,
+ nearestNeighbors(pressureGradient.centering(),
+ fo, disspoke, true)[0]);
+ // The column number is the pressure gradient corresponding to the
+ // "positive" side of the face.
+
+ // FIXME: The lhs (double) and rhs (vector field) do not match.
+
+ A[faceIndex][columnNu] =
+ directedPermeability(nearestNeighbors(directedPermeability.centering(),
+ fo, disspoke, true)[0]);
+ A[faceIndex+size][columnNu] =
+ faceDistance(nearestNeighbors(faceDistance.centering(),
+ fo, disspoke, true)[0]);
+ rhs[faceIndex+size] -=
+ pressure(nearestNeighbors(pressure.centering(),
+ fo, disspoke, true)[0]);
+
+ fo = disFluxPoints[faceIndex][1];
+ columnNu =
+ findIndex(gradients,
+ nearestNeighbors(pressureGradient.centering(),
+ fo, disspoke, true)[0]);
+ // The column number is the pressure gradient corresponding to the
+ // "positive" side of the face.
+
+ A[faceIndex][columnNu] =
+ -directedPermeability(nearestNeighbors(directedPermeability.centering(),
+ fo, disspoke, true)[0]);
+ A[faceIndex+size][columnNu] =
+ -faceDistance(nearestNeighbors(faceDistance.centering(),
+ fo, disspoke, true)[0]);
+ rhs[faceIndex+size] -=
+ -pressure(nearestNeighbors(pressure.centering(),
+ fo, disspoke, true)[0]);
+ }
! // Solve for the pressure gradients.
!
! TNT::LU_solve(A, ipiv, rhs);
!
! // Now, rhs has the pressure gradients.
! for (int faceIndex = size-1; faceIndex >= 0; --faceIndex)
! // FIXME: Is this type of assignment supported by the current code base?
! pressureGradient(gradients[faceIndex], subcell) = rhs[faceIndex];
!
!
! // Compute the spoke fluxes.
spokeFlux =
! dot(replicate(pressureGradient,
! nearestNeighbors(pressureGradient.centering(),
! spokeFlux.centering(),
! true),
! spokeFlux.centering()),
! replicate(directedPermeability,
! nearestNeighbors(directedPermeability.centering(),
! spokeFlux.centering(),
! true),
! spokeFlux.centering()));
// Sum the spoke fluxes into a cell flux.
*************** int main(int argc, char *argv[])
*** 217,224 ****
totalFlux =
sum(spokeFlux.mesh().normals().signedMagnitude() *
! sum(spokeFlux, nearestNeighbors(spoke, disFace)),
! nearestNeighbors(disFace, cell));
/* TERMINATION */
--- 346,364 ----
totalFlux =
sum(spokeFlux.mesh().normals().signedMagnitude() *
! // FIXME: This is not yet implemented. We want a
! // FIXME: data-parallel sum. This is we want a function
! // FIXME: Field_t sum(/* input */ Field_t,
! // FIXME: std::vector<FieldOffsetList>, /*output */
! // FIXME: Centering). The vector's length == the output
! // centering's length. The function works by using the input
! // field with each FieldOffsetList to form one value in the
! // output field.
! sum(spokeFlux,
! nearestNeighbors(spokeFlux.centering(), disFace),
! disFace),
! nearestNeighbors(disFace, totalFlux.centering()),
! totalFlux.centering());
/* TERMINATION */
Index: makefile
===================================================================
RCS file: /home/pooma/Repository/r2/examples/NewField/StatigraphicFlow/Attic/makefile,v
retrieving revision 1.1.2.1
diff -c -p -r1.1.2.1 makefile
*** makefile 2001/08/04 02:19:06 1.1.2.1
--- makefile 2001/08/16 20:26:07
*************** $(ODIR)/StatigraphicFlow: $(ODIR)/Statig
*** 44,49 ****
--- 44,51 ----
include $(SHARED_ROOT)/tail.mk
+ RULE_INCLUDES += -I$(PROJECT_ROOT)/examples/NewField/StatigraphicFlow # permit inclusion of the TNT linear algebra library
+
# ACL:rcsinfo
# ----------------------------------------------------------------------
# $RCSfile: makefile,v $ $Author: oldham $
Index: tnt/A10x10.dat
===================================================================
RCS file: A10x10.dat
diff -N A10x10.dat
*** /dev/null Tue May 5 14:32:27 1998
--- A10x10.dat Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,12 ----
+ 10 10
+ 0.218959 0.5297 0.526929 0.328234 0.0726859 0.166507 0.493977 0.464446 0.629543 0.845982
+ 0.0470446 0.671149 0.0919649 0.632639 0.631635 0.486517 0.266145 0.94098 0.736225 0.412081
+ 0.678865 0.00769819 0.653919 0.75641 0.884707 0.897656 0.0907329 0.050084 0.725412 0.841511
+ 0.679296 0.383416 0.415999 0.991037 0.27271 0.909208 0.947764 0.761514 0.999458 0.269317
+ 0.934693 0.0668422 0.701191 0.365339 0.436411 0.0605643 0.0737491 0.770205 0.888572 0.415395
+ 0.383502 0.417486 0.910321 0.247039 0.766495 0.904653 0.500707 0.827817 0.233195 0.537304
+ 0.519416 0.686773 0.762198 0.98255 0.477732 0.504523 0.384142 0.125365 0.306322 0.467917
+ 0.830965 0.588977 0.262453 0.72266 0.237774 0.516292 0.277082 0.0158677 0.351015 0.287212
+ 0.0345721 0.930436 0.0474645 0.753356 0.274907 0.319033 0.913817 0.688455 0.513274 0.178328
+ 0.0534616 0.846167 0.736082 0.651519 0.359265 0.986642 0.529747 0.868247 0.591114 0.15372
+
Index: tnt/AB10x10.dat
===================================================================
RCS file: AB10x10.dat
diff -N AB10x10.dat
*** /dev/null Tue May 5 14:32:27 1998
--- AB10x10.dat Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,24 ----
+ 10 10
+ 0.218959 0.5297 0.526929 0.328234 0.0726859 0.166507 0.493977 0.464446 0.629543 0.845982
+ 0.0470446 0.671149 0.0919649 0.632639 0.631635 0.486517 0.266145 0.94098 0.736225 0.412081
+ 0.678865 0.00769819 0.653919 0.75641 0.884707 0.897656 0.0907329 0.050084 0.725412 0.841511
+ 0.679296 0.383416 0.415999 0.991037 0.27271 0.909208 0.947764 0.761514 0.999458 0.269317
+ 0.934693 0.0668422 0.701191 0.365339 0.436411 0.0605643 0.0737491 0.770205 0.888572 0.415395
+ 0.383502 0.417486 0.910321 0.247039 0.766495 0.904653 0.500707 0.827817 0.233195 0.537304
+ 0.519416 0.686773 0.762198 0.98255 0.477732 0.504523 0.384142 0.125365 0.306322 0.467917
+ 0.830965 0.588977 0.262453 0.72266 0.237774 0.516292 0.277082 0.0158677 0.351015 0.287212
+ 0.0345721 0.930436 0.0474645 0.753356 0.274907 0.319033 0.913817 0.688455 0.513274 0.178328
+ 0.0534616 0.846167 0.736082 0.651519 0.359265 0.986642 0.529747 0.868247 0.591114 0.15372
+
+ 10 10
+ 0.571655 0.84204 0.70982 0.387725 0.408767 0.629269 0.901673 0.365339 0.215248 0.462245
+ 0.802406 0.159768 0.937897 0.499741 0.14182 0.126712 0.426497 0.253057 0.679592 0.951367
+ 0.0330538 0.212752 0.239911 0.147533 0.564899 0.651254 0.142021 0.135109 0.908922 0.632739
+ 0.53445 0.71471 0.180896 0.587187 0.252126 0.621634 0.947487 0.783153 0.250126 0.43933
+ 0.49848 0.130427 0.31754 0.845576 0.488515 0.803073 0.410313 0.455307 0.86086 0.824697
+ 0.955361 0.0909903 0.886991 0.590109 0.464031 0.247842 0.131189 0.349524 0.471262 0.688981
+ 0.748293 0.274588 0.652059 0.955409 0.961095 0.476432 0.885648 0.4523 0.505956 0.702207
+ 0.554584 0.0029996 0.150335 0.556146 0.126031 0.389314 0.0921736 0.808945 0.600394 0.987145
+ 0.890737 0.414293 0.681346 0.148152 0.199757 0.20325 0.162199 0.931674 0.817561 0.954415
+ 0.624849 0.0268763 0.385815 0.983305 0.31925 0.0283752 0.0710636 0.651646 0.755844 0.85127
+
Index: tnt/Makefile
===================================================================
RCS file: Makefile
diff -N Makefile
*** /dev/null Tue May 5 14:32:27 1998
--- Makefile Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,29 ----
+ ### Oldham, Jeffrey D.
+ ### 1999 Nov 22
+ ### programming
+ ###
+ ### Compile C++ Program
+
+ WCFLAGS= -pedantic -Wall -W -Wstrict-prototypes -Wpointer-arith -Wbad-function-cast -Wcast-align -Wconversion -Wnested-externs -Wundef -Winline
+ CFLAGS= -g $(WCFLAGS) -I.. # -DTIME -DSTATS
+ CC= ${HOME}/gcc-install/gcc1/bin/gcc
+ CXXFLAGS= -DUSE_LIBGXX_INLINES $(CFLAGS)
+ CXX= ${HOME}/gcc-install/gcc1/bin/g++
+ LDFLAGS= #-lm
+
+ CCSOURCES = main.cc
+
+ OBJECTS = $(CCSOURCES:%.cc=%.o)
+
+
+ %: %.o
+ $(CXX) $(CXXFLAGS) $^ $(LIBS) -o $@ $(LDFLAGS)
+
+ clean:
+ rm -f *.o
+
+ header-dependencies:
+ gcc -MM $(CFLAGS) $(CCSOURCES) $(CSOURCES)
+
+ ## ADD header file dependencies
+ ## Create them using "gmake -k header-dependencies".
Index: tnt/SPD10x10.dat
===================================================================
RCS file: SPD10x10.dat
diff -N SPD10x10.dat
*** /dev/null Tue May 5 14:32:27 1998
--- SPD10x10.dat Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,11 ----
+ 10 10
+ 2.3186576e+00 2.1294424e+00 2.1960039e+00 2.7464350e+00 2.0762250e+00 2.3053886e+00 2.1570514e+00 1.5808710e+00 2.0909577e+00 2.4191432e+00
+ 2.1294424e+00 3.1651186e+00 2.5234253e+00 3.3847286e+00 2.2593036e+00 2.7678316e+00 2.3628136e+00 1.7825304e+00 2.7783113e+00 3.2137616e+00
+ 2.1960039e+00 2.5234253e+00 4.2942805e+00 3.6189729e+00 2.8497810e+00 3.2440901e+00 3.1321684e+00 2.4829696e+00 1.8009041e+00 2.8701790e+00
+ 2.7464350e+00 3.3847286e+00 3.6189729e+00 5.2143007e+00 3.1447822e+00 3.5583692e+00 3.3876913e+00 2.8527866e+00 3.4629890e+00 4.1031277e+00
+ 2.0762250e+00 2.2593036e+00 2.8497810e+00 3.1447822e+00 3.2581495e+00 2.6091394e+00 2.2553001e+00 1.8630124e+00 1.6701147e+00 2.3741290e+00
+ 2.3053886e+00 2.7678316e+00 3.2440901e+00 3.5583692e+00 2.6091394e+00 3.8960566e+00 2.8640541e+00 2.0193737e+00 2.3733266e+00 3.5771674e+00
+ 2.1570514e+00 2.3628136e+00 3.1321684e+00 3.3876913e+00 2.2553001e+00 2.8640541e+00 3.2466334e+00 2.4706183e+00 2.4036480e+00 3.0448446e+00
+ 1.5808710e+00 1.7825304e+00 2.4829696e+00 2.8527866e+00 1.8630124e+00 2.0193737e+00 2.4706183e+00 2.2343394e+00 1.8592024e+00 2.2138309e+00
+ 2.0909577e+00 2.7783113e+00 1.8009041e+00 3.4629890e+00 1.6701147e+00 2.3733266e+00 2.4036480e+00 1.8592024e+00 3.2183447e+00 3.1411090e+00
+ 2.4191432e+00 3.2137616e+00 2.8701790e+00 4.1031277e+00 2.3741290e+00 3.5771674e+00 3.0448446e+00 2.2138309e+00 3.1411090e+00 4.1952140e+00
Index: tnt/cholesky.h
===================================================================
RCS file: cholesky.h
diff -N cholesky.h
*** /dev/null Tue May 5 14:32:27 1998
--- cholesky.h Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,96 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ #ifndef CHOLESKY_H
+ #define CHOLESKY_H
+
+ #include <cmath>
+
+ // index method
+
+ namespace TNT
+ {
+
+
+ //
+ // Only upper part of A is used. Cholesky factor is returned in
+ // lower part of L. Returns 0 if successful, 1 otherwise.
+ //
+ template <class SPDMatrix, class SymmMatrix>
+ int Cholesky_upper_factorization(SPDMatrix &A, SymmMatrix &L)
+ {
+ Subscript M = A.dim(1);
+ Subscript N = A.dim(2);
+
+ assert(M == N); // make sure A is square
+
+ // readjust size of L, if necessary
+
+ if (M != L.dim(1) || N != L.dim(2))
+ L = SymmMatrix(N,N);
+
+ Subscript i,j,k;
+
+
+ typename SPDMatrix::element_type dot=0;
+
+
+ for (j=1; j<=N; j++) // form column j of L
+ {
+ dot= 0;
+
+ for (i=1; i<j; i++) // for k= 1 TO j-1
+ dot = dot + L(j,i)*L(j,i);
+
+ L(j,j) = A(j,j) - dot;
+
+ for (i=j+1; i<=N; i++)
+ {
+ dot = 0;
+ for (k=1; k<j; k++)
+ dot = dot + L(i,k)*L(j,k);
+ L(i,j) = A(j,i) - dot;
+ }
+
+ if (L(j,j) <= 0.0) return 1;
+
+ L(j,j) = sqrt( L(j,j) );
+
+ for (i=j+1; i<=N; i++)
+ L(i,j) = L(i,j) / L(j,j);
+
+ }
+
+ return 0;
+ }
+
+
+
+
+ }
+ // namespace TNT
+
+ #endif
+ // CHOLESKY_H
Index: tnt/cmat.h
===================================================================
RCS file: cmat.h
diff -N cmat.h
*** /dev/null Tue May 5 14:32:27 1998
--- cmat.h Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,605 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ // C compatible matrix: row-oriented, 0-based [i][j] and 1-based (i,j) indexing
+ //
+
+ #ifndef CMAT_H
+ #define CMAT_H
+
+ #include "tnt/subscript.h"
+ #include "tnt/vec.h"
+ #include <cstdlib>
+ #include <cassert>
+ #include <iostream>
+ #include <sstream>
+ #ifdef TNT_USE_REGIONS
+ #include "tnt/region2d.h"
+ #endif
+
+ namespace TNT
+ {
+
+ template <class T>
+ class Matrix
+ {
+
+
+ public:
+
+ typedef Subscript size_type;
+ typedef T value_type;
+ typedef T element_type;
+ typedef T* pointer;
+ typedef T* iterator;
+ typedef T& reference;
+ typedef const T* const_iterator;
+ typedef const T& const_reference;
+
+ Subscript lbound() const { return 1;}
+
+ protected:
+ Subscript m_;
+ Subscript n_;
+ Subscript mn_; // total size
+ T* v_;
+ T** row_;
+ T* vm1_ ; // these point to the same data, but are 1-based
+ T** rowm1_;
+
+ // internal helper function to create the array
+ // of row pointers
+
+ void initialize(Subscript M, Subscript N)
+ {
+ mn_ = M*N;
+ m_ = M;
+ n_ = N;
+
+ v_ = new T[mn_];
+ row_ = new T*[M];
+ rowm1_ = new T*[M];
+
+ assert(v_ != NULL);
+ assert(row_ != NULL);
+ assert(rowm1_ != NULL);
+
+ T* p = v_;
+ vm1_ = v_ - 1;
+ for (Subscript i=0; i<M; i++)
+ {
+ row_[i] = p;
+ rowm1_[i] = p-1;
+ p += N ;
+
+ }
+
+ rowm1_ -- ; // compensate for 1-based offset
+ }
+
+ void copy(const T* v)
+ {
+ Subscript N = m_ * n_;
+ Subscript i;
+
+ #ifdef TNT_UNROLL_LOOPS
+ Subscript Nmod4 = N & 3;
+ Subscript N4 = N - Nmod4;
+
+ for (i=0; i<N4; i+=4)
+ {
+ v_[i] = v[i];
+ v_[i+1] = v[i+1];
+ v_[i+2] = v[i+2];
+ v_[i+3] = v[i+3];
+ }
+
+ for (i=N4; i< N; i++)
+ v_[i] = v[i];
+ #else
+
+ for (i=0; i< N; i++)
+ v_[i] = v[i];
+ #endif
+ }
+
+ void set(const T& val)
+ {
+ Subscript N = m_ * n_;
+ Subscript i;
+
+ #ifdef TNT_UNROLL_LOOPS
+ Subscript Nmod4 = N & 3;
+ Subscript N4 = N - Nmod4;
+
+ for (i=0; i<N4; i+=4)
+ {
+ v_[i] = val;
+ v_[i+1] = val;
+ v_[i+2] = val;
+ v_[i+3] = val;
+ }
+
+ for (i=N4; i< N; i++)
+ v_[i] = val;
+ #else
+
+ for (i=0; i< N; i++)
+ v_[i] = val;
+
+ #endif
+ }
+
+
+
+ void destroy()
+ {
+ /* do nothing, if no memory has been previously allocated */
+ if (v_ == NULL) return ;
+
+ /* if we are here, then matrix was previously allocated */
+ if (v_ != NULL) delete [] (v_);
+ if (row_ != NULL) delete [] (row_);
+
+ /* return rowm1_ back to original value */
+ rowm1_ ++;
+ if (rowm1_ != NULL ) delete [] (rowm1_);
+ }
+
+
+ public:
+
+ operator T**(){ return row_; }
+ operator T**() const { return row_; }
+
+
+ Subscript size() const { return mn_; }
+
+ // constructors
+
+ Matrix() : m_(0), n_(0), mn_(0), v_(0), row_(0), vm1_(0), rowm1_(0) {};
+
+ Matrix(const Matrix<T> &A)
+ {
+ initialize(A.m_, A.n_);
+ copy(A.v_);
+ }
+
+ Matrix(Subscript M, Subscript N, const T& value = T())
+ {
+ initialize(M,N);
+ set(value);
+ }
+
+ Matrix(Subscript M, Subscript N, const T* v)
+ {
+ initialize(M,N);
+ copy(v);
+ }
+
+ Matrix(Subscript M, Subscript N, const char *s)
+ {
+ initialize(M,N);
+ std::istringstream ins(s);
+
+ Subscript i, j;
+
+ for (i=0; i<M; i++)
+ for (j=0; j<N; j++)
+ ins >> row_[i][j];
+ }
+
+ // destructor
+ //
+ ~Matrix()
+ {
+ destroy();
+ }
+
+
+ // reallocating
+ //
+ Matrix<T>& newsize(Subscript M, Subscript N)
+ {
+ if (num_rows() == M && num_cols() == N)
+ return *this;
+
+ destroy();
+ initialize(M,N);
+
+ return *this;
+ }
+
+
+
+
+ // assignments
+ //
+ Matrix<T>& operator=(const Matrix<T> &A)
+ {
+ if (v_ == A.v_)
+ return *this;
+
+ if (m_ == A.m_ && n_ == A.n_) // no need to re-alloc
+ copy(A.v_);
+
+ else
+ {
+ destroy();
+ initialize(A.m_, A.n_);
+ copy(A.v_);
+ }
+
+ return *this;
+ }
+
+ Matrix<T>& operator=(const T& scalar)
+ {
+ set(scalar);
+ return *this;
+ }
+
+
+ Subscript dim(Subscript d) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( d >= 1);
+ assert( d <= 2);
+ #endif
+ return (d==1) ? m_ : ((d==2) ? n_ : 0);
+ }
+
+ Subscript num_rows() const { return m_; }
+ Subscript num_cols() const { return n_; }
+
+
+
+
+ inline T* operator[](Subscript i)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(0<=i);
+ assert(i < m_) ;
+ #endif
+ return row_[i];
+ }
+
+ inline const T* operator[](Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(0<=i);
+ assert(i < m_) ;
+ #endif
+ return row_[i];
+ }
+
+ inline reference operator()(Subscript i)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i <= mn_) ;
+ #endif
+ return vm1_[i];
+ }
+
+ inline const_reference operator()(Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i <= mn_) ;
+ #endif
+ return vm1_[i];
+ }
+
+
+
+ inline reference operator()(Subscript i, Subscript j)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i <= m_) ;
+ assert(1<=j);
+ assert(j <= n_);
+ #endif
+ return rowm1_[i][j];
+ }
+
+
+
+ inline const_reference operator() (Subscript i, Subscript j) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i <= m_) ;
+ assert(1<=j);
+ assert(j <= n_);
+ #endif
+ return rowm1_[i][j];
+ }
+
+
+
+ #ifdef TNT_USE_REGIONS
+
+ typedef Region2D<Matrix<T> > Region;
+
+
+ Region operator()(const Index1D &I, const Index1D &J)
+ {
+ return Region(*this, I,J);
+ }
+
+
+ typedef const_Region2D< Matrix<T> > const_Region;
+ const_Region operator()(const Index1D &I, const Index1D &J) const
+ {
+ return const_Region(*this, I,J);
+ }
+
+ #endif
+
+
+ };
+
+
+ /* *************************** I/O ********************************/
+
+ template <class T>
+ std::ostream& operator<<(std::ostream &s, const Matrix<T> &A)
+ {
+ Subscript M=A.num_rows();
+ Subscript N=A.num_cols();
+
+ s << M << " " << N << "\n";
+
+ for (Subscript i=0; i<M; i++)
+ {
+ for (Subscript j=0; j<N; j++)
+ {
+ s << A[i][j] << " ";
+ }
+ s << "\n";
+ }
+
+
+ return s;
+ }
+
+ template <class T>
+ std::istream& operator>>(std::istream &s, Matrix<T> &A)
+ {
+
+ Subscript M, N;
+
+ s >> M >> N;
+
+ if ( !(M == A.num_rows() && N == A.num_cols() ))
+ {
+ A.newsize(M,N);
+ }
+
+
+ for (Subscript i=0; i<M; i++)
+ for (Subscript j=0; j<N; j++)
+ {
+ s >> A[i][j];
+ }
+
+
+ return s;
+ }
+
+ // *******************[ basic matrix algorithms ]***************************
+
+
+ template <class T>
+ Matrix<T> operator+(const Matrix<T> &A,
+ const Matrix<T> &B)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ assert(M==B.num_rows());
+ assert(N==B.num_cols());
+
+ Matrix<T> tmp(M,N);
+ Subscript i,j;
+
+ for (i=0; i<M; i++)
+ for (j=0; j<N; j++)
+ tmp[i][j] = A[i][j] + B[i][j];
+
+ return tmp;
+ }
+
+ template <class T>
+ Matrix<T> operator-(const Matrix<T> &A,
+ const Matrix<T> &B)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ assert(M==B.num_rows());
+ assert(N==B.num_cols());
+
+ Matrix<T> tmp(M,N);
+ Subscript i,j;
+
+ for (i=0; i<M; i++)
+ for (j=0; j<N; j++)
+ tmp[i][j] = A[i][j] - B[i][j];
+
+ return tmp;
+ }
+
+ template <class T>
+ Matrix<T> mult_element(const Matrix<T> &A,
+ const Matrix<T> &B)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ assert(M==B.num_rows());
+ assert(N==B.num_cols());
+
+ Matrix<T> tmp(M,N);
+ Subscript i,j;
+
+ for (i=0; i<M; i++)
+ for (j=0; j<N; j++)
+ tmp[i][j] = A[i][j] * B[i][j];
+
+ return tmp;
+ }
+
+
+ template <class T>
+ Matrix<T> transpose(const Matrix<T> &A)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ Matrix<T> S(N,M);
+ Subscript i, j;
+
+ for (i=0; i<M; i++)
+ for (j=0; j<N; j++)
+ S[j][i] = A[i][j];
+
+ return S;
+ }
+
+
+
+ template <class T>
+ inline Matrix<T> matmult(const Matrix<T> &A,
+ const Matrix<T> &B)
+ {
+
+ #ifdef TNT_BOUNDS_CHECK
+ assert(A.num_cols() == B.num_rows());
+ #endif
+
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+ Subscript K = B.num_cols();
+
+ Matrix<T> tmp(M,K);
+ T sum;
+
+ for (Subscript i=0; i<M; i++)
+ for (Subscript k=0; k<K; k++)
+ {
+ sum = 0;
+ for (Subscript j=0; j<N; j++)
+ sum = sum + A[i][j] * B[j][k];
+
+ tmp[i][k] = sum;
+ }
+
+ return tmp;
+ }
+
+ template <class T>
+ inline Matrix<T> operator*(const Matrix<T> &A,
+ const Matrix<T> &B)
+ {
+ return matmult(A,B);
+ }
+
+ template <class T>
+ inline int matmult(Matrix<T>& C, const Matrix<T> &A,
+ const Matrix<T> &B)
+ {
+
+ assert(A.num_cols() == B.num_rows());
+
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+ Subscript K = B.num_cols();
+
+ C.newsize(M,K);
+
+ T sum;
+
+ const T* row_i;
+ const T* col_k;
+
+ for (Subscript i=0; i<M; i++)
+ for (Subscript k=0; k<K; k++)
+ {
+ row_i = &(A[i][0]);
+ col_k = &(B[0][k]);
+ sum = 0;
+ for (Subscript j=0; j<N; j++)
+ {
+ sum += *row_i * *col_k;
+ row_i++;
+ col_k += K;
+ }
+ C[i][k] = sum;
+ }
+
+ return 0;
+ }
+
+
+ template <class T>
+ Vector<T> matmult(const Matrix<T> &A, const Vector<T> &x)
+ {
+
+ #ifdef TNT_BOUNDS_CHECK
+ assert(A.num_cols() == x.dim());
+ #endif
+
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ Vector<T> tmp(M);
+ T sum;
+
+ for (Subscript i=0; i<M; i++)
+ {
+ sum = 0;
+ const T* rowi = A[i];
+ for (Subscript j=0; j<N; j++)
+ sum = sum + rowi[j] * x[j];
+
+ tmp[i] = sum;
+ }
+
+ return tmp;
+ }
+
+ template <class T>
+ inline Vector<T> operator*(const Matrix<T> &A, const Vector<T> &x)
+ {
+ return matmult(A,x);
+ }
+
+ } // namespace TNT
+
+ #endif
+ // CMAT_H
Index: tnt/fchol.cc
===================================================================
RCS file: fchol.cc
diff -N fchol.cc
*** /dev/null Tue May 5 14:32:27 1998
--- fchol.cc Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,60 ----
+
+
+
+ // Test Cholesky module
+
+ #include <iostream>
+
+ #include "tnt/tnt.h"
+ #include "tnt/vec.h"
+ #include "tnt/fmat.h"
+ #include "tnt/cholesky.h"
+ #include "tnt/trisolve.h"
+ #include "tnt/transv.h" /* transpose views */
+
+ using namespace std;
+ using namespace TNT;
+
+ int main()
+ {
+ Fortran_Matrix<double> A;
+
+ cin >> A; /* A should be symmetric positive definite */
+
+ Subscript N = A.num_rows();
+ assert(N == A.num_cols());
+
+ Vector<double> b(N, 1.0); // b= [1,1,1,...]
+ Fortran_Matrix<double> L(N, N);
+
+
+ cout << "A: " << A << endl;
+
+ if (Cholesky_upper_factorization(A, L) !=0)
+ {
+ cout << "Cholesky did not work." << endl;
+ exit(1);
+ }
+
+
+ cout << L << endl;
+
+ // solve Ax =b, as L*L'x =b
+ //
+ // let y=L'x, then
+ //
+ //
+ // solve L y = b;
+ // solve L'x = y;
+
+ Vector<double> y = Lower_triangular_solve(L, b);
+ TNT::Transpose_View<TNT::Fortran_Matrix<double> > foo(L);
+ Vector<double> x1=
+ Upper_triangular_solve(Transpose_View<Fortran_Matrix<double> >(L), y);
+ Vector<double> x= Upper_triangular_solve(foo, y);
+
+ cout << "x: " << x << endl;
+ cout << "Residual A*x-b: " << A*x-b << endl;
+
+ return 0;
+ }
Index: tnt/fcscmat.h
===================================================================
RCS file: fcscmat.h
diff -N fcscmat.h
*** /dev/null Tue May 5 14:32:27 1998
--- fcscmat.h Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,165 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ // Templated compressed sparse column matrix (Fortran conventions).
+ // uses 1-based offsets in storing row indices.
+ // Used primarily to interface with Fortran sparse matrix libaries.
+ // (CANNOT BE USED AS AN STL CONTAINER.)
+
+
+ #ifndef FCSCMAT_H
+ #define FCSCMAT_H
+
+ #include <iostream>
+ #include <cassert>
+ #include "tnt/tnt.h"
+ #include "tnt/vec.h"
+
+ using namespace std;
+
+ namespace TNT
+ {
+
+ template <class T>
+ class Fortran_Sparse_Col_Matrix
+ {
+
+ protected:
+
+ Vector<T> val_; // data values (nz_ elements)
+ Vector<Subscript> rowind_; // row_ind (nz_ elements)
+ Vector<Subscript> colptr_; // col_ptr (n_+1 elements)
+
+ int nz_; // number of nonzeros
+ Subscript m_; // global dimensions
+ Subscript n_;
+
+ public:
+
+
+ Fortran_Sparse_Col_Matrix(void);
+ Fortran_Sparse_Col_Matrix(const Fortran_Sparse_Col_Matrix<T> &S)
+ : val_(S.val_), rowind_(S.rowind_), colptr_(S.colptr_), nz_(S.nz_),
+ m_(S.m_), n_(S.n_) {};
+ Fortran_Sparse_Col_Matrix(Subscript M, Subscript N,
+ Subscript nz, const T *val, const Subscript *r,
+ const Subscript *c) : val_(nz, val), rowind_(nz, r),
+ colptr_(N+1, c), nz_(nz), m_(M), n_(N) {};
+
+ Fortran_Sparse_Col_Matrix(Subscript M, Subscript N,
+ Subscript nz, char *val, char *r,
+ char *c) : val_(nz, val), rowind_(nz, r),
+ colptr_(N+1, c), nz_(nz), m_(M), n_(N) {};
+
+ Fortran_Sparse_Col_Matrix(Subscript M, Subscript N,
+ Subscript nz, const T *val, Subscript *r, Subscript *c)
+ : val_(nz, val), rowind_(nz, r), colptr_(N+1, c), nz_(nz),
+ m_(M), n_(N) {};
+
+ ~Fortran_Sparse_Col_Matrix() {};
+
+
+ T & val(Subscript i) { return val_(i); }
+ const T & val(Subscript i) const { return val_(i); }
+
+ Subscript & row_ind(Subscript i) { return rowind_(i); }
+ const Subscript & row_ind(Subscript i) const { return rowind_(i); }
+
+ Subscript col_ptr(Subscript i) { return colptr_(i);}
+ const Subscript col_ptr(Subscript i) const { return colptr_(i);}
+
+
+ Subscript num_cols() const { return m_;}
+ Subscript num_rows() const { return n_; }
+
+ Subscript dim(Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( 1 <= i );
+ assert( i <= 2 );
+ #endif
+ if (i==1) return m_;
+ else if (i==2) return m_;
+ else return 0;
+ }
+
+ Subscript num_nonzeros() const {return nz_;};
+ Subscript lbound() const {return 1;}
+
+
+
+ Fortran_Sparse_Col_Matrix& operator=(const
+ Fortran_Sparse_Col_Matrix &C)
+ {
+ val_ = C.val_;
+ rowind_ = C.rowind_;
+ colptr_ = C.colptr_;
+ nz_ = C.nz_;
+ m_ = C.m_;
+ n_ = C.n_;
+
+ return *this;
+ }
+
+ Fortran_Sparse_Col_Matrix& newsize(Subscript M, Subscript N,
+ Subscript nz)
+ {
+ val_.newsize(nz);
+ rowind_.newsize(nz);
+ colptr_.newsize(N+1);
+ return *this;
+ }
+ };
+
+ template <class T>
+ ostream& operator<<(ostream &s, const Fortran_Sparse_Col_Matrix<T> &A)
+ {
+ Subscript M=A.num_rows();
+ Subscript N=A.num_cols();
+
+ s << M << " " << N << " " << A.num_nonzeros() << endl;
+
+
+ for (Subscript k=1; k<=N; k++)
+ {
+ Subscript start = A.col_ptr(k);
+ Subscript end = A.col_ptr(k+1);
+
+ for (Subscript i= start; i<end; i++)
+ {
+ s << A.row_ind(i) << " " << k << " " << A.val(i) << endl;
+ }
+ }
+
+ return s;
+ }
+
+
+
+ } // namespace TNT
+
+ #endif
+ /* FCSCMAT_H */
+
Index: tnt/fmat.h
===================================================================
RCS file: fmat.h
diff -N fmat.h
*** /dev/null Tue May 5 14:32:27 1998
--- fmat.h Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,577 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ // Fortran-compatible matrix: column oriented, 1-based (i,j) indexing
+
+ #ifndef FMAT_H
+ #define FMAT_H
+
+ #include "tnt/subscript.h"
+ #include "tnt/vec.h"
+ #include <cstdlib>
+ #include <cassert>
+ #include <iostream>
+ #include <sstream>
+ #ifdef TNT_USE_REGIONS
+ #include "tnt/region2d.h"
+ #endif
+
+ // simple 1-based, column oriented Matrix class
+
+ namespace TNT
+ {
+
+ template <class T>
+ class Fortran_Matrix
+ {
+
+
+ public:
+
+ typedef T value_type;
+ typedef T element_type;
+ typedef T* pointer;
+ typedef T* iterator;
+ typedef T& reference;
+ typedef const T* const_iterator;
+ typedef const T& const_reference;
+
+ Subscript lbound() const { return 1;}
+
+ protected:
+ T* v_; // these are adjusted to simulate 1-offset
+ Subscript m_;
+ Subscript n_;
+ T** col_; // these are adjusted to simulate 1-offset
+
+ // internal helper function to create the array
+ // of row pointers
+
+ void initialize(Subscript M, Subscript N)
+ {
+ // adjust col_[] pointers so that they are 1-offset:
+ // col_[j][i] is really col_[j-1][i-1];
+ //
+ // v_[] is the internal contiguous array, it is still 0-offset
+ //
+ v_ = new T[M*N];
+ col_ = new T*[N];
+
+ assert(v_ != NULL);
+ assert(col_ != NULL);
+
+
+ m_ = M;
+ n_ = N;
+ T* p = v_ - 1;
+ for (Subscript i=0; i<N; i++)
+ {
+ col_[i] = p;
+ p += M ;
+
+ }
+ col_ --;
+ }
+
+ void copy(const T* v)
+ {
+ Subscript N = m_ * n_;
+ Subscript i;
+
+ #ifdef TNT_UNROLL_LOOPS
+ Subscript Nmod4 = N & 3;
+ Subscript N4 = N - Nmod4;
+
+ for (i=0; i<N4; i+=4)
+ {
+ v_[i] = v[i];
+ v_[i+1] = v[i+1];
+ v_[i+2] = v[i+2];
+ v_[i+3] = v[i+3];
+ }
+
+ for (i=N4; i< N; i++)
+ v_[i] = v[i];
+ #else
+
+ for (i=0; i< N; i++)
+ v_[i] = v[i];
+ #endif
+ }
+
+ void set(const T& val)
+ {
+ Subscript N = m_ * n_;
+ Subscript i;
+
+ #ifdef TNT_UNROLL_LOOPS
+ Subscript Nmod4 = N & 3;
+ Subscript N4 = N - Nmod4;
+
+ for (i=0; i<N4; i+=4)
+ {
+ v_[i] = val;
+ v_[i+1] = val;
+ v_[i+2] = val;
+ v_[i+3] = val;
+ }
+
+ for (i=N4; i< N; i++)
+ v_[i] = val;
+ #else
+
+ for (i=0; i< N; i++)
+ v_[i] = val;
+
+ #endif
+ }
+
+
+
+ void destroy()
+ {
+ /* do nothing, if no memory has been previously allocated */
+ if (v_ == NULL) return ;
+
+ /* if we are here, then matrix was previously allocated */
+ delete [] (v_);
+ col_ ++; // changed back to 0-offset
+ delete [] (col_);
+ }
+
+
+ public:
+
+ T* begin() { return v_; }
+ const T* begin() const { return v_;}
+
+ T* end() { return v_ + m_*n_; }
+ const T* end() const { return v_ + m_*n_; }
+
+
+ // constructors
+
+ Fortran_Matrix() : v_(0), m_(0), n_(0), col_(0) {};
+ Fortran_Matrix(const Fortran_Matrix<T> &A)
+ {
+ initialize(A.m_, A.n_);
+ copy(A.v_);
+ }
+
+ Fortran_Matrix(Subscript M, Subscript N, const T& value = T())
+ {
+ initialize(M,N);
+ set(value);
+ }
+
+ Fortran_Matrix(Subscript M, Subscript N, const T* v)
+ {
+ initialize(M,N);
+ copy(v);
+ }
+
+
+ Fortran_Matrix(Subscript M, Subscript N, char *s)
+ {
+ initialize(M,N);
+ std::istringstream ins(s);
+
+ Subscript i, j;
+
+ for (i=1; i<=M; i++)
+ for (j=1; j<=N; j++)
+ ins >> (*this)(i,j);
+ }
+
+ // destructor
+ ~Fortran_Matrix()
+ {
+ destroy();
+ }
+
+
+ // assignments
+ //
+ Fortran_Matrix<T>& operator=(const Fortran_Matrix<T> &A)
+ {
+ if (v_ == A.v_)
+ return *this;
+
+ if (m_ == A.m_ && n_ == A.n_) // no need to re-alloc
+ copy(A.v_);
+
+ else
+ {
+ destroy();
+ initialize(A.m_, A.n_);
+ copy(A.v_);
+ }
+
+ return *this;
+ }
+
+ Fortran_Matrix<T>& operator=(const T& scalar)
+ {
+ set(scalar);
+ return *this;
+ }
+
+
+ Subscript dim(Subscript d) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( d >= 1);
+ assert( d <= 2);
+ #endif
+ return (d==1) ? m_ : ((d==2) ? n_ : 0);
+ }
+
+ Subscript num_rows() const { return m_; }
+ Subscript num_cols() const { return n_; }
+
+ Fortran_Matrix<T>& newsize(Subscript M, Subscript N)
+ {
+ if (num_rows() == M && num_cols() == N)
+ return *this;
+
+ destroy();
+ initialize(M,N);
+
+ return *this;
+ }
+
+
+
+ // 1-based element access
+ //
+ inline reference operator()(Subscript i, Subscript j)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i <= m_) ;
+ assert(1<=j);
+ assert(j <= n_);
+ #endif
+ return col_[j][i];
+ }
+
+ inline const_reference operator() (Subscript i, Subscript j) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i <= m_) ;
+ assert(1<=j);
+ assert(j <= n_);
+ #endif
+ return col_[j][i];
+ }
+
+
+ #ifdef TNT_USE_REGIONS
+
+ typedef Region2D<Fortran_Matrix<T> > Region;
+ typedef const_Region2D< Fortran_Matrix<T> > const_Region;
+
+ Region operator()(const Index1D &I, const Index1D &J)
+ {
+ return Region(*this, I,J);
+ }
+
+ const_Region operator()(const Index1D &I, const Index1D &J) const
+ {
+ return const_Region(*this, I,J);
+ }
+
+ #endif
+
+
+ };
+
+
+ /* *************************** I/O ********************************/
+
+ template <class T>
+ std::ostream& operator<<(std::ostream &s, const Fortran_Matrix<T> &A)
+ {
+ Subscript M=A.num_rows();
+ Subscript N=A.num_cols();
+
+ s << M << " " << N << "\n";
+
+ for (Subscript i=1; i<=M; i++)
+ {
+ for (Subscript j=1; j<=N; j++)
+ {
+ s << A(i,j) << " ";
+ }
+ s << "\n";
+ }
+
+
+ return s;
+ }
+
+ template <class T>
+ std::istream& operator>>(std::istream &s, Fortran_Matrix<T> &A)
+ {
+
+ Subscript M, N;
+
+ s >> M >> N;
+
+ if ( !(M == A.num_rows() && N == A.num_cols()))
+ {
+ A.newsize(M,N);
+ }
+
+
+ for (Subscript i=1; i<=M; i++)
+ for (Subscript j=1; j<=N; j++)
+ {
+ s >> A(i,j);
+ }
+
+
+ return s;
+ }
+
+ // *******************[ basic matrix algorithms ]***************************
+
+
+ template <class T>
+ Fortran_Matrix<T> operator+(const Fortran_Matrix<T> &A,
+ const Fortran_Matrix<T> &B)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ assert(M==B.num_rows());
+ assert(N==B.num_cols());
+
+ Fortran_Matrix<T> tmp(M,N);
+ Subscript i,j;
+
+ for (i=1; i<=M; i++)
+ for (j=1; j<=N; j++)
+ tmp(i,j) = A(i,j) + B(i,j);
+
+ return tmp;
+ }
+
+ template <class T>
+ Fortran_Matrix<T> operator-(const Fortran_Matrix<T> &A,
+ const Fortran_Matrix<T> &B)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ assert(M==B.num_rows());
+ assert(N==B.num_cols());
+
+ Fortran_Matrix<T> tmp(M,N);
+ Subscript i,j;
+
+ for (i=1; i<=M; i++)
+ for (j=1; j<=N; j++)
+ tmp(i,j) = A(i,j) - B(i,j);
+
+ return tmp;
+ }
+
+ // element-wise multiplication (use matmult() below for matrix
+ // multiplication in the linear algebra sense.)
+ //
+ //
+ template <class T>
+ Fortran_Matrix<T> mult_element(const Fortran_Matrix<T> &A,
+ const Fortran_Matrix<T> &B)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ assert(M==B.num_rows());
+ assert(N==B.num_cols());
+
+ Fortran_Matrix<T> tmp(M,N);
+ Subscript i,j;
+
+ for (i=1; i<=M; i++)
+ for (j=1; j<=N; j++)
+ tmp(i,j) = A(i,j) * B(i,j);
+
+ return tmp;
+ }
+
+
+ template <class T>
+ Fortran_Matrix<T> transpose(const Fortran_Matrix<T> &A)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ Fortran_Matrix<T> S(N,M);
+ Subscript i, j;
+
+ for (i=1; i<=M; i++)
+ for (j=1; j<=N; j++)
+ S(j,i) = A(i,j);
+
+ return S;
+ }
+
+
+
+ template <class T>
+ inline Fortran_Matrix<T> matmult(const Fortran_Matrix<T> &A,
+ const Fortran_Matrix<T> &B)
+ {
+
+ #ifdef TNT_BOUNDS_CHECK
+ assert(A.num_cols() == B.num_rows());
+ #endif
+
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+ Subscript K = B.num_cols();
+
+ Fortran_Matrix<T> tmp(M,K);
+ T sum;
+
+ for (Subscript i=1; i<=M; i++)
+ for (Subscript k=1; k<=K; k++)
+ {
+ sum = 0;
+ for (Subscript j=1; j<=N; j++)
+ sum = sum + A(i,j) * B(j,k);
+
+ tmp(i,k) = sum;
+ }
+
+ return tmp;
+ }
+
+ template <class T>
+ inline Fortran_Matrix<T> operator*(const Fortran_Matrix<T> &A,
+ const Fortran_Matrix<T> &B)
+ {
+ return matmult(A,B);
+ }
+
+ template <class T>
+ inline int matmult(Fortran_Matrix<T>& C, const Fortran_Matrix<T> &A,
+ const Fortran_Matrix<T> &B)
+ {
+
+ assert(A.num_cols() == B.num_rows());
+
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+ Subscript K = B.num_cols();
+
+ C.newsize(M,K); // adjust shape of C, if necessary
+
+
+ T sum;
+
+ const T* row_i;
+ const T* col_k;
+
+ for (Subscript i=1; i<=M; i++)
+ {
+ for (Subscript k=1; k<=K; k++)
+ {
+ row_i = &A(i,1);
+ col_k = &B(1,k);
+ sum = 0;
+ for (Subscript j=1; j<=N; j++)
+ {
+ sum += *row_i * *col_k;
+ row_i += M;
+ col_k ++;
+ }
+
+ C(i,k) = sum;
+ }
+
+ }
+
+ return 0;
+ }
+
+
+ template <class T>
+ Vector<T> matmult(const Fortran_Matrix<T> &A, const Vector<T> &x)
+ {
+
+ #ifdef TNT_BOUNDS_CHECK
+ assert(A.num_cols() == x.dim());
+ #endif
+
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ Vector<T> tmp(M);
+ T sum;
+
+ for (Subscript i=1; i<=M; i++)
+ {
+ sum = 0;
+ for (Subscript j=1; j<=N; j++)
+ sum = sum + A(i,j) * x(j);
+
+ tmp(i) = sum;
+ }
+
+ return tmp;
+ }
+
+ template <class T>
+ inline Vector<T> operator*(const Fortran_Matrix<T> &A, const Vector<T> &x)
+ {
+ return matmult(A,x);
+ }
+
+ template <class T>
+ inline Fortran_Matrix<T> operator*(const Fortran_Matrix<T> &A, const T &x)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ Subscript MN = M*N;
+
+ Fortran_Matrix<T> res(M,N);
+ const T* a = A.begin();
+ T* t = res.begin();
+ T* tend = res.end();
+
+ for (t=res.begin(); t < tend; t++, a++)
+ *t = *a * x;
+
+ return res;
+ }
+
+ } // namespace TNT
+ #endif
+ // FMAT_H
Index: tnt/fortran.h
===================================================================
RCS file: fortran.h
diff -N fortran.h
*** /dev/null Tue May 5 14:32:27 1998
--- fortran.h Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,67 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ // Header file to define C/Fortran conventions (Platform specific)
+
+ #ifndef FORTRAN_H
+ #define FORTRAN_H
+
+ // help map between C/C++ data types and Fortran types
+
+ typedef int Fortran_integer;
+ typedef float Fortran_float;
+ typedef double Fortran_double;
+
+
+ typedef Fortran_double *fda_; // (in/out) double precision array
+ typedef const Fortran_double *cfda_; // (in) double precsion array
+
+ typedef Fortran_double *fd_; // (in/out) single double precision
+ typedef const Fortran_double *cfd_; // (in) single double precision
+
+ typedef Fortran_float *ffa_; // (in/out) float precision array
+ typedef const Fortran_float *cffa_; // (in) float precsion array
+
+ typedef Fortran_float *ff_; // (in/out) single float precision
+ typedef const Fortran_float *cff_; // (in) single float precision
+
+ typedef Fortran_integer *fia_; // (in/out) single integer array
+ typedef const Fortran_integer *cfia_; // (in) single integer array
+
+ typedef Fortran_integer *fi_; // (in/out) single integer
+ typedef const Fortran_integer *cfi_; // (in) single integer
+
+ typedef char *fch_; // (in/out) single character
+ typedef char *cfch_; // (in) single character
+
+
+
+ #ifndef TNT_SUBSCRIPT_TYPE
+ #define TNT_SUBSCRIPT_TYPE TNT::Fortran_integer
+ #endif
+
+
+ #endif
+ // FORTRAN_H
Index: tnt/fspvec.h
===================================================================
RCS file: fspvec.h
diff -N fspvec.h
*** /dev/null Tue May 5 14:32:27 1998
--- fspvec.h Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,168 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+ // Templated sparse vector (Fortran conventions).
+ // Used primarily to interface with Fortran sparse matrix libaries.
+ // (CANNOT BE USED AS AN STL CONTAINER.)
+
+ #ifndef FSPVEC_H
+ #define FSPVEC_H
+
+ #include "tnt/tnt.h"
+ #include "tnt/vec.h"
+ #include <cstdlib>
+ #include <cassert>
+ #include <iostream>
+ #include <sstream>
+
+ using namespace std;
+
+ namespace TNT
+ {
+
+ template <class T>
+ class Fortran_Sparse_Vector
+ {
+
+
+ public:
+
+ typedef Subscript size_type;
+ typedef T value_type;
+ typedef T element_type;
+ typedef T* pointer;
+ typedef T* iterator;
+ typedef T& reference;
+ typedef const T* const_iterator;
+ typedef const T& const_reference;
+
+ Subscript lbound() const { return 1;}
+
+ protected:
+ Vector<T> val_;
+ Vector<Subscript> index_;
+ Subscript dim_; // prescribed dimension
+
+
+ public:
+
+ // size and shape information
+
+ Subscript dim() const { return dim_; }
+ Subscript num_nonzeros() const { return val_.dim(); }
+
+ // access
+
+ T& val(Subscript i) { return val_(i); }
+ const T& val(Subscript i) const { return val_(i); }
+
+ Subscript &index(Subscript i) { return index_(i); }
+ const Subscript &index(Subscript i) const { return index_(i); }
+
+ // constructors
+
+ Fortran_Sparse_Vector() : val_(), index_(), dim_(0) {};
+ Fortran_Sparse_Vector(Subscript N, Subscript nz) : val_(nz),
+ index_(nz), dim_(N) {};
+ Fortran_Sparse_Vector(Subscript N, Subscript nz, const T *values,
+ const Subscript *indices): val_(nz, values), index_(nz, indices),
+ dim_(N) {}
+
+ Fortran_Sparse_Vector(const Fortran_Sparse_Vector<T> &S):
+ val_(S.val_), index_(S.index_), dim_(S.dim_) {}
+
+ // initialize from string, e.g.
+ //
+ // Fortran_Sparse_Vector<T> A(N, 2, "1.0 2.1", "1 3");
+ //
+ Fortran_Sparse_Vector(Subscript N, Subscript nz, char *v,
+ char *ind) : val_(nz, v), index_(nz, ind), dim_(N) {}
+
+ // assignments
+
+ Fortran_Sparse_Vector<T> & newsize(Subscript N, Subscript nz)
+ {
+ val_.newsize(nz);
+ index_.newsize(nz);
+ dim_ = N;
+ return *this;
+ }
+
+ Fortran_Sparse_Vector<T> & operator=( const Fortran_Sparse_Vector<T> &A)
+ {
+ val_ = A.val_;
+ index_ = A.index_;
+ dim_ = A.dim_;
+
+ return *this;
+ }
+
+ // methods
+
+
+
+ };
+
+
+ /* *************************** I/O ********************************/
+
+ template <class T>
+ ostream& operator<<(ostream &s, const Fortran_Sparse_Vector<T> &A)
+ {
+ // output format is : N nz val1 ind1 val2 ind2 ...
+ Subscript nz=A.num_nonzeros();
+
+ s << A.dim() << " " << nz << endl;
+
+ for (Subscript i=1; i<=nz; i++)
+ s << A.val(i) << " " << A.index(i) << endl;
+ s << endl;
+
+ return s;
+ }
+
+
+ template <class T>
+ istream& operator>>(istream &s, Fortran_Sparse_Vector<T> &A)
+ {
+ // output format is : N nz val1 ind1 val2 ind2 ...
+
+ Subscript N;
+ Subscript nz;
+
+ s >> N >> nz;
+
+ A.newsize(N, nz);
+
+ for (Subscript i=1; i<=nz; i++)
+ s >> A.val(i) >> A.index(i);
+
+
+ return s;
+ }
+
+ } // namespace TNT
+
+ #endif
+ // FSPVEC_H
Index: tnt/index.h
===================================================================
RCS file: index.h
diff -N index.h
*** /dev/null Tue May 5 14:32:27 1998
--- index.h Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,83 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ // Vector/Matrix/Array Index Module
+
+ #ifndef INDEX_H
+ #define INDEX_H
+
+ #include "tnt/subscript.h"
+
+ namespace TNT
+ {
+
+ class Index1D
+ {
+ Subscript lbound_;
+ Subscript ubound_;
+
+ public:
+
+ Subscript lbound() const { return lbound_; }
+ Subscript ubound() const { return ubound_; }
+
+ Index1D(const Index1D &D) : lbound_(D.lbound_), ubound_(D.ubound_) {}
+ Index1D(Subscript i1, Subscript i2) : lbound_(i1), ubound_(i2) {}
+
+ Index1D & operator=(const Index1D &D)
+ {
+ lbound_ = D.lbound_;
+ ubound_ = D.ubound_;
+ return *this;
+ }
+
+ };
+
+ inline Index1D operator+(const Index1D &D, Subscript i)
+ {
+ return Index1D(i+D.lbound(), i+D.ubound());
+ }
+
+ inline Index1D operator+(Subscript i, const Index1D &D)
+ {
+ return Index1D(i+D.lbound(), i+D.ubound());
+ }
+
+
+
+ inline Index1D operator-(Index1D &D, Subscript i)
+ {
+ return Index1D(D.lbound()-i, D.ubound()-i);
+ }
+
+ inline Index1D operator-(Subscript i, Index1D &D)
+ {
+ return Index1D(i-D.lbound(), i-D.ubound());
+ }
+
+ } // namespace TNT
+
+ #endif
+
Index: tnt/lapack.h
===================================================================
RCS file: lapack.h
diff -N lapack.h
*** /dev/null Tue May 5 14:32:27 1998
--- lapack.h Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,193 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ // Header file for Fortran Lapack
+
+ #ifndef LAPACK_H
+ #define LAPACK_H
+
+ // This file incomplete and included here to only demonstrate the
+ // basic framework for linking with the Fortran Lapack routines.
+
+ #include "tnt/fortran.h"
+ #include "tnt/vec.h"
+ #include "tnt/fmat.h"
+
+
+ #define F77_DGESV dgesv_
+ #define F77_DGELS dgels_
+ #define F77_DSYEV dsyev_
+ #define F77_DGEEV dgeev_
+
+ extern "C"
+ {
+
+ // linear equations (general) using LU factorizaiton
+ //
+ void F77_DGESV(cfi_ N, cfi_ nrhs, fda_ A, cfi_ lda,
+ fia_ ipiv, fda_ b, cfi_ ldb, fi_ info);
+
+ // solve linear least squares using QR or LU factorization
+ //
+ void F77_DGELS(cfch_ trans, cfi_ M,
+ cfi_ N, cfi_ nrhs, fda_ A, cfi_ lda, fda_ B, cfi_ ldb, fda_ work,
+ cfi_ lwork, fi_ info);
+
+ // solve symmetric eigenvalues
+ //
+ void F77_DSYEV( cfch_ jobz, cfch_ uplo, cfi_ N, fda_ A, cfi_ lda,
+ fda_ W, fda_ work, cfi_ lwork, fi_ info);
+
+ // solve unsymmetric eigenvalues
+ //
+ void F77_DGEEV(cfch_ jobvl, cfch_ jobvr, cfi_ N, fda_ A, cfi_ lda,
+ fda_ wr, fda_ wi, fda_ vl, cfi_ ldvl, fda_ vr,
+ cfi_ ldvr, fda_ work, cfi_ lwork, fi_ info);
+
+ }
+
+ // solve linear equations using LU factorization
+
+ using namespace TNT;
+
+ Vector<double> Lapack_LU_linear_solve(const Fortran_Matrix<double> &A,
+ const Vector<double> &b)
+ {
+ const Fortran_integer one=1;
+ Subscript M=A.num_rows();
+ Subscript N=A.num_cols();
+
+ Fortran_Matrix<double> Tmp(A);
+ Vector<double> x(b);
+ Vector<Fortran_integer> index(M);
+ Fortran_integer info = 0;
+
+ F77_DGESV(&N, &one, &Tmp(1,1), &M, &index(1), &x(1), &M, &info);
+
+ if (info != 0) return Vector<double>(0);
+ else
+ return x;
+ }
+
+ // solve linear least squares problem using QR factorization
+ //
+ Vector<double> Lapack_LLS_QR_linear_solve(const Fortran_Matrix<double> &A,
+ const Vector<double> &b)
+ {
+ const Fortran_integer one=1;
+ Subscript M=A.num_rows();
+ Subscript N=A.num_cols();
+
+ Fortran_Matrix<double> Tmp(A);
+ Vector<double> x(b);
+ Fortran_integer info = 0;
+
+ char transp = 'N';
+ Fortran_integer lwork = 5 * (M+N); // temporary work space
+ Vector<double> work(lwork);
+
+ F77_DGELS(&transp, &M, &N, &one, &Tmp(1,1), &M, &x(1), &M, &work(1),
+ &lwork, &info);
+
+ if (info != 0) return Vector<double>(0);
+ else
+ return x;
+ }
+
+ // *********************** Eigenvalue problems *******************
+
+ // solve symmetric eigenvalue problem (eigenvalues only)
+ //
+ Vector<double> Upper_symmetric_eigenvalue_solve(const Fortran_Matrix<double> &A)
+ {
+ char jobz = 'N';
+ char uplo = 'U';
+ Subscript N = A.num_rows();
+
+ assert(N == A.num_cols());
+
+ Vector<double> eigvals(N);
+ Fortran_integer worksize = 3*N;
+ Fortran_integer info = 0;
+ Vector<double> work(worksize);
+ Fortran_Matrix<double> Tmp = A;
+
+ F77_DSYEV(&jobz, &uplo, &N, &Tmp(1,1), &N, eigvals.begin(), work.begin(),
+ &worksize, &info);
+
+ if (info != 0) return Vector<double>();
+ else
+ return eigvals;
+ }
+
+
+ // solve unsymmetric eigenvalue problems
+ //
+ int eigenvalue_solve(const Fortran_Matrix<double> &A,
+ Vector<double> &wr, Vector<double> &wi)
+ {
+ char jobvl = 'N';
+ char jobvr = 'N';
+
+ Fortran_integer N = A.num_rows();
+
+
+ assert(N == A.num_cols());
+
+ if (N<1) return 1;
+
+ Fortran_Matrix<double> vl(1,N); /* should be NxN ? **** */
+ Fortran_Matrix<double> vr(1,N);
+ Fortran_integer one = 1;
+
+ Fortran_integer worksize = 5*N;
+ Fortran_integer info = 0;
+ Vector<double> work(worksize, 0.0);
+ Fortran_Matrix<double> Tmp = A;
+
+ wr.newsize(N);
+ wi.newsize(N);
+
+ // void F77_DGEEV(cfch_ jobvl, cfch_ jobvr, cfi_ N, fda_ A, cfi_ lda,
+ // fda_ wr, fda_ wi, fda_ vl, cfi_ ldvl, fda_ vr,
+ // cfi_ ldvr, fda_ work, cfi_ lwork, fi_ info);
+
+ F77_DGEEV(&jobvl, &jobvr, &N, &Tmp(1,1), &N, &(wr(1)),
+ &(wi(1)), &(vl(1,1)), &one, &(vr(1,1)), &one,
+ &(work(1)), &worksize, &info);
+
+ return (info==0 ? 0: 1);
+ }
+
+
+
+
+
+ #endif
+ // LAPACK_H
+
+
+
+
Index: tnt/lu-C.cc
===================================================================
RCS file: lu-C.cc
diff -N lu-C.cc
*** /dev/null Tue May 5 14:32:27 1998
--- lu-C.cc Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,57 ----
+ // Solve a linear system using LU factorization.
+ //
+ // Usage: a.out < matrix.dat
+ //
+ // where matrix.dat is an ASCII file consisting of the
+ // matrix size (M,N) followed by its values. For example,
+ //
+ // 3 3
+ // 8.1 1.2 4.3
+ // 1.3 4.3 2.9
+ // 0.4 1.3 6.1
+
+ #include <iostream>
+
+ #include "tnt/tnt.h"
+ #include "tnt/vec.h"
+ #include "tnt/cmat.h"
+ #include "tnt/lu.h"
+
+ using namespace std;
+ using namespace TNT;
+
+ int main()
+ {
+ Matrix<double> A;
+
+ cin >> A;
+
+ Subscript N = A.dim(1);
+ assert(N == A.dim(2));
+
+ Vector<double> b(N, 1.0); // b= [1,1,1,...]
+ Vector<Subscript> index(N);
+
+ cout << "Original Matrix A: " << A << endl;
+
+ Matrix<double> T(A);
+ if (LU_factor(T, index) !=0)
+ {
+ cout << "LU_factor() failed." << endl;
+ exit(1);
+ }
+
+ Vector<double> x(b);
+ if (LU_solve(T, index, x) != 0)
+ {
+ cout << "LU_Solve() failed." << endl;
+ exit(1);
+ }
+ cout << "Solution x for Ax=b, where b=[1,1,...] " <<endl;
+ cout << " x: " << x << endl;
+
+ cout << "A*x should be the vector [1,1,...] " <<endl;
+ cout << "residual [A*x - b]: " << matmult(A, x) - b << endl;
+
+ return 0;
+ }
Index: tnt/lu.cc
===================================================================
RCS file: lu.cc
diff -N lu.cc
*** /dev/null Tue May 5 14:32:27 1998
--- lu.cc Thu Aug 16 13:02:55 2001
***************
*** 0 ****
--- 1,64 ----
+
+
+ // Solve a linear system using LU factorization.
+ //
+ // Usage: a.out < matrix.dat
+ //
+ // where matrix.dat is an ASCII file consisting of the
+ // matrix size (M,N) followed by its values. For example,
+ //
+ // 3 3
+ // 8.1 1.2 4.3
+ // 1.3 4.3 2.9
+ // 0.4 1.3 6.1
+ //
+
+
+ #include <iostream>
+
+ #include "tnt/tnt.h"
+ #include "tnt/vec.h"
+ #include "tnt/fmat.h"
+ #include "tnt/lu.h"
+
+ using namespace std;
+ using namespace TNT;
+
+ int main()
+ {
+ Fortran_Matrix<double> A;
+
+ cin >> A;
+
+
+ Subscript N = A.dim(1);
+ assert(N == A.dim(2));
+
+ Vector<double> b(N, 1.0); // b= [1,1,1,...]
+ Vector<Subscript> index(N);
+
+
+
+ cout << "Original Matrix A: " << A << endl;
+
+ Fortran_Matrix<double> T(A);
+ if (LU_factor(T, index) !=0)
+ {
+ cout << "LU_factor() failed." << endl;
+ exit(1);
+ }
+
+ Vector<double> x(b);
+ if (LU_solve(T, index, x) != 0)
+ {
+ cout << "LU_Solve() failed." << endl;
+ exit(1);
+ }
+ cout << "Solution x for Ax=b, where b=[1,1,...] " <<endl;
+ cout << " x: " << x << endl;
+
+ cout << "A*x should be the vector [1,1,...] " <<endl;
+ cout << "residual [A*x - b]: " << matmult(A, x) - b << endl;
+
+ return 0;
+ }
Index: tnt/lu.h
===================================================================
RCS file: lu.h
diff -N lu.h
*** /dev/null Tue May 5 14:32:27 1998
--- lu.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,204 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ #ifndef LU_H
+ #define LU_H
+
+ // Solve system of linear equations Ax = b.
+ //
+ // Typical usage:
+ //
+ // Matrix<double> A;
+ // Vector<Subscript> ipiv;
+ // Vector<double> b;
+ //
+ // 1) LU_factor(A,ipiv);
+ // 2) LU_solve(A,ipiv,b);
+ //
+ // Now b has the solution x. Note that both A and b
+ // are overwritten. If these values need to be preserved,
+ // one can make temporary copies, as in
+ //
+ // O) Matrix<double> T = A;
+ // 1) LU_factor(T,ipiv);
+ // 1a) Vector<double> x=b;
+ // 2) LU_solve(T,ipiv,x);
+ //
+ // See details below.
+ //
+
+
+ // for fabs()
+ //
+ #include <cmath>
+
+ // right-looking LU factorization algorithm (unblocked)
+ //
+ // Factors matrix A into lower and upper triangular matrices
+ // (L and U respectively) in solving the linear equation Ax=b.
+ //
+ //
+ // Args:
+ //
+ // A (input/output) Matrix(1:n, 1:n) In input, matrix to be
+ // factored. On output, overwritten with lower and
+ // upper triangular factors.
+ //
+ // indx (output) Vector(1:n) Pivot vector. Describes how
+ // the rows of A were reordered to increase
+ // numerical stability.
+ //
+ // Return value:
+ //
+ // int (0 if successful, 1 otherwise)
+ //
+ //
+
+
+ namespace TNT
+ {
+
+ template <class MaTRiX, class VecToRSubscript>
+ int LU_factor( MaTRiX &A, VecToRSubscript &indx)
+ {
+ assert(A.lbound() == 1); // currently for 1-offset
+ assert(indx.lbound() == 1); // vectors and matrices
+
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ if (M == 0 || N==0) return 0;
+ if (indx.dim() != M)
+ indx.newsize(M);
+
+ Subscript i=0,j=0,k=0;
+ Subscript jp=0;
+
+ typename MaTRiX::element_type t;
+
+ Subscript minMN = (M < N ? M : N) ; // min(M,N);
+
+ for (j=1; j<= minMN; j++)
+ {
+
+ // find pivot in column j and test for singularity.
+
+ jp = j;
+ t = fabs(A(j,j));
+ for (i=j+1; i<=M; i++)
+ if ( fabs(A(i,j)) > t)
+ {
+ jp = i;
+ t = fabs(A(i,j));
+ }
+
+ indx(j) = jp;
+
+ // jp now has the index of maximum element
+ // of column j, below the diagonal
+
+ if ( A(jp,j) == 0 )
+ return 1; // factorization failed because of zero pivot
+
+
+ if (jp != j) // swap rows j and jp
+ for (k=1; k<=N; k++)
+ {
+ t = A(j,k);
+ A(j,k) = A(jp,k);
+ A(jp,k) =t;
+ }
+
+ if (j<M) // compute elements j+1:M of jth column
+ {
+ // note A(j,j), was A(jp,p) previously which was
+ // guarranteed not to be zero (Label #1)
+ //
+ typename MaTRiX::element_type recp = 1.0 / A(j,j);
+
+ for (k=j+1; k<=M; k++)
+ A(k,j) *= recp;
+ }
+
+
+ if (j < minMN)
+ {
+ // rank-1 update to trailing submatrix: E = E - x*y;
+ //
+ // E is the region A(j+1:M, j+1:N)
+ // x is the column vector A(j+1:M,j)
+ // y is row vector A(j,j+1:N)
+
+ Subscript ii,jj;
+
+ for (ii=j+1; ii<=M; ii++)
+ for (jj=j+1; jj<=N; jj++)
+ A(ii,jj) -= A(ii,j)*A(j,jj);
+ }
+ }
+
+ return 0;
+ }
+
+
+
+
+ template <class MaTRiX, class VecToR, class VecToRSubscripts>
+ int LU_solve(const MaTRiX &A, const VecToRSubscripts &indx, VecToR &b)
+ {
+ assert(A.lbound() == 1); // currently for 1-offset
+ assert(indx.lbound() == 1); // vectors and matrices
+ assert(b.lbound() == 1);
+
+ Subscript i,ii=0,ip,j;
+ Subscript n = b.dim();
+ typename MaTRiX::element_type sum = 0.0;
+
+ for (i=1;i<=n;i++)
+ {
+ ip=indx(i);
+ sum=b(ip);
+ b(ip)=b(i);
+ if (ii)
+ for (j=ii;j<=i-1;j++)
+ sum -= A(i,j)*b(j);
+ else if (sum) ii=i;
+ b(i)=sum;
+ }
+ for (i=n;i>=1;i--)
+ {
+ sum=b(i);
+ for (j=i+1;j<=n;j++)
+ sum -= A(i,j)*b(j);
+ b(i)=sum/A(i,i);
+ }
+
+ return 0;
+ }
+
+ } // namespace TNT
+
+ #endif
+ // LU_H
Index: tnt/matrix.dat
===================================================================
RCS file: matrix.dat
diff -N matrix.dat
*** /dev/null Tue May 5 14:32:27 1998
--- matrix.dat Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,4 ----
+ 3 3
+ 8.1 1.2 4.3
+ 1.3 4.3 2.9
+ 0.4 1.3 6.1
Index: tnt/qr.cc
===================================================================
RCS file: qr.cc
diff -N qr.cc
*** /dev/null Tue May 5 14:32:27 1998
--- qr.cc Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,66 ----
+
+
+
+ // Solve a linear system using QR factorization.
+ //
+ // Usage: a.out < matrix.dat
+ //
+ // where matrix.dat is an ASCII file consisting of the
+ // matrix size (M,N) followed by its values. For example,
+ //
+ // 3 2
+ // 8.1 1.2 4.3
+ // 1.3 4.3 2.9
+ // 0.4 1.3 6.1
+ //
+ //
+
+ #include <iostream>
+
+ #include "tnt/tnt.h"
+ #include "tnt/vec.h"
+ #include "tnt/fmat.h"
+ #include "tnt/qr.h"
+
+ using namespace std;
+ using namespace TNT;
+
+ int main()
+
+ {
+ Fortran_Matrix<double> A;
+
+ cin >> A;
+
+ Subscript N = A.num_rows();
+ assert(N == A.num_cols());
+
+ Vector<double> b(N, 1.0); // b= [1,1,1,...]
+ Vector<double> C(N), D(N);
+
+
+ cout << "A: " << A << endl;
+
+ Fortran_Matrix<double> T(A);
+
+ if (QR_factor(T, C, D) !=0)
+ {
+ cout << "QR failed." << endl;
+ cout << " returned: \n" << T << endl;
+ exit(1);
+ }
+
+ Vector<double> x(b);
+ if (QR_solve(T, C, D, x) == 1)
+ {
+ cout << "QR_Solve did not work." << endl;
+ exit(1);
+ }
+
+ cout << "Solution x for Ax=b, where b=[1,1,...] " <<endl;
+ cout << " x: " << x << endl;
+
+ cout << "residual [A*x - b]: " << A*x - b << endl;
+
+ return 0;
+ }
Index: tnt/qr.h
===================================================================
RCS file: qr.h
diff -N qr.h
*** /dev/null Tue May 5 14:32:27 1998
--- qr.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,229 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+ #ifndef QR_H
+ #define QR_H
+
+ // Classical QR factorization example, based on Stewart[1973].
+ //
+ //
+ // This algorithm computes the factorization of a matrix A
+ // into a product of an orthognal matrix (Q) and an upper triangular
+ // matrix (R), such that QR = A.
+ //
+ // Parameters:
+ //
+ // A (in): Matrix(1:N, 1:N)
+ //
+ // Q (output): Matrix(1:N, 1:N), collection of Householder
+ // column vectors Q1, Q2, ... QN
+ //
+ // R (output): upper triangular Matrix(1:N, 1:N)
+ //
+ // Returns:
+ //
+ // 0 if successful, 1 if A is detected to be singular
+ //
+
+
+ #include <cmath> //for sqrt() & fabs()
+ #include "tnt/tntmath.h" // for sign()
+
+ // Classical QR factorization, based on Stewart[1973].
+ //
+ //
+ // This algorithm computes the factorization of a matrix A
+ // into a product of an orthognal matrix (Q) and an upper triangular
+ // matrix (R), such that QR = A.
+ //
+ // Parameters:
+ //
+ // A (in/out): On input, A is square, Matrix(1:N, 1:N), that represents
+ // the matrix to be factored.
+ //
+ // On output, Q and R is encoded in the same Matrix(1:N,1:N)
+ // in the following manner:
+ //
+ // R is contained in the upper triangular section of A,
+ // except that R's main diagonal is in D. The lower
+ // triangular section of A represents Q, where each
+ // column j is the vector Qj = I - uj*uj'/pi_j.
+ //
+ // C (output): vector of Pi[j]
+ // D (output): main diagonal of R, i.e. D(i) is R(i,i)
+ //
+ // Returns:
+ //
+ // 0 if successful, 1 if A is detected to be singular
+ //
+
+ namespace TNT
+ {
+
+ template <class MaTRiX, class Vector>
+ int QR_factor(MaTRiX &A, Vector& C, Vector &D)
+ {
+ assert(A.lbound() == 1); // ensure these are all
+ assert(C.lbound() == 1); // 1-based arrays and vectors
+ assert(D.lbound() == 1);
+
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ assert(M == N); // make sure A is square
+
+ Subscript i,j,k;
+ typename MaTRiX::element_type eta, sigma, sum;
+
+ // adjust the shape of C and D, if needed...
+
+ if (N != C.size()) C.newsize(N);
+ if (N != D.size()) D.newsize(N);
+
+ for (k=1; k<N; k++)
+ {
+ // eta = max |M(i,k)|, for k <= i <= n
+ //
+ eta = 0;
+ for (i=k; i<=N; i++)
+ {
+ double absA = fabs(A(i,k));
+ eta = ( absA > eta ? absA : eta );
+ }
+
+ if (eta == 0) // matrix is singular
+ {
+ cerr << "QR: k=" << k << "\n";
+ return 1;
+ }
+
+ // form Qk and premiltiply M by it
+ //
+ for(i=k; i<=N; i++)
+ A(i,k) = A(i,k) / eta;
+
+ sum = 0;
+ for (i=k; i<=N; i++)
+ sum = sum + A(i,k)*A(i,k);
+ sigma = sign(A(k,k)) * sqrt(sum);
+
+
+ A(k,k) = A(k,k) + sigma;
+ C(k) = sigma * A(k,k);
+ D(k) = -eta * sigma;
+
+ for (j=k+1; j<=N; j++)
+ {
+ sum = 0;
+ for (i=k; i<=N; i++)
+ sum = sum + A(i,k)*A(i,j);
+ sum = sum / C(k);
+
+ for (i=k; i<=N; i++)
+ A(i,j) = A(i,j) - sum * A(i,k);
+ }
+
+ D(N) = A(N,N);
+ }
+
+ return 0;
+ }
+
+ // modified form of upper triangular solve, except that the main diagonal
+ // of R (upper portion of A) is in D.
+ //
+ template <class MaTRiX, class Vector>
+ int R_solve(const MaTRiX &A, /*const*/ Vector &D, Vector &b)
+ {
+ assert(A.lbound() == 1); // ensure these are all
+ assert(D.lbound() == 1); // 1-based arrays and vectors
+ assert(b.lbound() == 1);
+
+ Subscript i,j;
+ Subscript N = A.num_rows();
+
+ assert(N == A.num_cols());
+ assert(N == D.dim());
+ assert(N == b.dim());
+
+ typename MaTRiX::element_type sum;
+
+ if (D(N) == 0)
+ return 1;
+
+ b(N) = b(N) /
+ D(N);
+
+ for (i=N-1; i>=1; i--)
+ {
+ if (D(i) == 0)
+ return 1;
+ sum = 0;
+ for (j=i+1; j<=N; j++)
+ sum = sum + A(i,j)*b(j);
+ b(i) = ( b(i) - sum ) /
+ D(i);
+ }
+
+ return 0;
+ }
+
+
+ template <class MaTRiX, class Vector>
+ int QR_solve(const MaTRiX &A, const Vector &c, /*const*/ Vector &d,
+ Vector &b)
+ {
+ assert(A.lbound() == 1); // ensure these are all
+ assert(c.lbound() == 1); // 1-based arrays and vectors
+ assert(d.lbound() == 1);
+
+ Subscript N=A.num_rows();
+
+ assert(N == A.num_cols());
+ assert(N == c.dim());
+ assert(N == d.dim());
+ assert(N == b.dim());
+
+ Subscript i,j;
+ typename MaTRiX::element_type sum, tau;
+
+ for (j=1; j<N; j++)
+ {
+ // form Q'*b
+ sum = 0;
+ for (i=j; i<=N; i++)
+ sum = sum + A(i,j)*b(i);
+ if (c(j) == 0)
+ return 1;
+ tau = sum / c(j);
+ for (i=j; i<=N; i++)
+ b(i) = b(i) - tau * A(i,j);
+ }
+ return R_solve(A, d, b); // solve Rx = Q'b
+ }
+
+ } // namespace TNT
+
+ #endif
+ // QR_H
Index: tnt/region1d.h
===================================================================
RCS file: region1d.h
diff -N region1d.h
*** /dev/null Tue May 5 14:32:27 1998
--- region1d.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,371 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+
+ #ifndef REGION1D_H
+ #define REGION1D_H
+
+
+ #include "tnt/subscript.h"
+ #include "tnt/index.h"
+ #include <iostream>
+ #include <cassert>
+
+ namespace TNT
+ {
+
+ template <class Array1D>
+ class const_Region1D;
+
+ template <class Array1D>
+ class Region1D
+ {
+ protected:
+
+ Array1D & A_;
+ Subscript offset_; // 0-based
+ Subscript dim_;
+
+ typedef typename Array1D::element_type T;
+
+ public:
+ const Array1D & array() const { return A_; }
+
+ Subscript offset() const { return offset_;}
+ Subscript dim() const { return dim_; }
+
+ Subscript offset(Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(i==TNT_BASE_OFFSET);
+ #endif
+ return offset_;
+ }
+
+ Subscript dim(Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(i== TNT_BASE_OFFSET);
+ #endif
+ return offset_;
+ }
+
+
+ Region1D(Array1D &A, Subscript i1, Subscript i2) : A_(A)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(TNT_BASE_OFFSET <= i1 );
+ assert(i2 <= A.dim() + (TNT_BASE_OFFSET-1));
+ assert(i1 <= i2);
+ #endif
+ offset_ = i1 - TNT_BASE_OFFSET;
+ dim_ = i2-i1 + 1;
+ }
+
+ Region1D(Array1D &A, const Index1D &I) : A_(A)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(TNT_BASE_OFFSET <=I.lbound());
+ assert(I.ubound() <= A.dim() + (TNT_BASE_OFFSET-1));
+ assert(I.lbound() <= I.ubound());
+ #endif
+ offset_ = I.lbound() - TNT_BASE_OFFSET;
+ dim_ = I.ubound() - I.lbound() + 1;
+ }
+
+ Region1D(Region1D<Array1D> &A, Subscript i1, Subscript i2) :
+ A_(A.A_)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(TNT_BASE_OFFSET <= i1 );
+ assert(i2 <= A.dim() + (TNT_BASE_OFFSET - 1));
+ assert(i1 <= i2);
+ #endif
+ // (old-offset) (new-offset)
+ //
+ offset_ = (i1 - TNT_BASE_OFFSET) + A.offset_;
+ dim_ = i2-i1 + 1;
+ }
+
+ Region1D<Array1D> operator()(Subscript i1, Subscript i2)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(TNT_BASE_OFFSET <= i1);
+ assert(i2 <= dim() + (TNT_BASE_OFFSET -1));
+ assert(i1 <= i2);
+ #endif
+ // offset_ is 0-based, so no need for
+ // ( - TNT_BASE_OFFSET)
+ //
+ return Region1D<Array1D>(A_, i1+offset_,
+ offset_ + i2);
+ }
+
+
+ Region1D<Array1D> operator()(const Index1D &I)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(TNT_BASE_OFFSET<=I.lbound());
+ assert(I.ubound() <= dim() + (TNT_BASE_OFFSET-1));
+ assert(I.lbound() <= I.ubound());
+ #endif
+ return Region1D<Array1D>(A_, I.lbound()+offset_,
+ offset_ + I.ubound());
+ }
+
+
+
+
+ T & operator()(Subscript i)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(TNT_BASE_OFFSET <= i);
+ assert(i <= dim() + (TNT_BASE_OFFSET-1));
+ #endif
+ return A_(i+offset_);
+ }
+
+ const T & operator() (Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(TNT_BASE_OFFSET <= i);
+ assert(i <= dim() + (TNT_BASE_OFFSET-1));
+ #endif
+ return A_(i+offset_);
+ }
+
+
+ Region1D<Array1D> & operator=(const Region1D<Array1D> &R)
+ {
+ // make sure both sides conform
+ assert(dim() == R.dim());
+
+ Subscript N = dim();
+ Subscript i;
+ Subscript istart = TNT_BASE_OFFSET;
+ Subscript iend = istart + N-1;
+
+ for (i=istart; i<=iend; i++)
+ (*this)(i) = R(i);
+
+ return *this;
+ }
+
+
+
+ Region1D<Array1D> & operator=(const const_Region1D<Array1D> &R)
+ {
+ // make sure both sides conform
+ assert(dim() == R.dim());
+
+ Subscript N = dim();
+ Subscript i;
+ Subscript istart = TNT_BASE_OFFSET;
+ Subscript iend = istart + N-1;
+
+ for (i=istart; i<=iend; i++)
+ (*this)(i) = R(i);
+
+ return *this;
+
+ }
+
+
+ Region1D<Array1D> & operator=(const T& t)
+ {
+ Subscript N=dim();
+ Subscript i;
+ Subscript istart = TNT_BASE_OFFSET;
+ Subscript iend = istart + N-1;
+
+ for (i=istart; i<= iend; i++)
+ (*this)(i) = t;
+
+ return *this;
+
+ }
+
+
+ Region1D<Array1D> & operator=(const Array1D &R)
+ {
+ // make sure both sides conform
+ Subscript N = dim();
+ assert(dim() == R.dim());
+
+ Subscript i;
+ Subscript istart = TNT_BASE_OFFSET;
+ Subscript iend = istart + N-1;
+
+ for (i=istart; i<=iend; i++)
+ (*this)(i) = R(i);
+
+ return *this;
+
+ }
+
+ };
+
+ template <class Array1D>
+ std::ostream& operator<<(std::ostream &s, Region1D<Array1D> &A)
+ {
+ Subscript N=A.dim();
+ Subscript istart = TNT_BASE_OFFSET;
+ Subscript iend = N - 1 + TNT_BASE_OFFSET;
+
+ for (Subscript i=istart; i<=iend; i++)
+ s << A(i) << endl;
+
+ return s;
+ }
+
+
+ /* --------- class const_Region1D ------------ */
+
+ template <class Array1D>
+ class const_Region1D
+ {
+ protected:
+
+ const Array1D & A_;
+ Subscript offset_; // 0-based
+ Subscript dim_;
+ typedef typename Array1D::element_type T;
+
+ public:
+ const Array1D & array() const { return A_; }
+
+ Subscript offset() const { return offset_;}
+ Subscript dim() const { return dim_; }
+
+ Subscript offset(Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(i==TNT_BASE_OFFSET);
+ #endif
+ return offset_;
+ }
+
+ Subscript dim(Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(i== TNT_BASE_OFFSET);
+ #endif
+ return offset_;
+ }
+
+
+ const_Region1D(const Array1D &A, Subscript i1, Subscript i2) : A_(A)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(TNT_BASE_OFFSET <= i1 );
+ assert(i2 <= A.dim() + (TNT_BASE_OFFSET-1));
+ assert(i1 <= i2);
+ #endif
+ offset_ = i1 - TNT_BASE_OFFSET;
+ dim_ = i2-i1 + 1;
+ }
+
+ const_Region1D(const Array1D &A, const Index1D &I) : A_(A)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(TNT_BASE_OFFSET <=I.lbound());
+ assert(I.ubound() <= A.dim() + (TNT_BASE_OFFSET-1));
+ assert(I.lbound() <= I.ubound());
+ #endif
+ offset_ = I.lbound() - TNT_BASE_OFFSET;
+ dim_ = I.ubound() - I.lbound() + 1;
+ }
+
+ const_Region1D(const_Region1D<Array1D> &A, Subscript i1, Subscript i2) :
+ A_(A.A_)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(TNT_BASE_OFFSET <= i1 );
+ assert(i2 <= A.dim() + (TNT_BASE_OFFSET - 1));
+ assert(i1 <= i2);
+ #endif
+ // (old-offset) (new-offset)
+ //
+ offset_ = (i1 - TNT_BASE_OFFSET) + A.offset_;
+ dim_ = i2-i1 + 1;
+ }
+
+ const_Region1D<Array1D> operator()(Subscript i1, Subscript i2)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(TNT_BASE_OFFSET <= i1);
+ assert(i2 <= dim() + (TNT_BASE_OFFSET -1));
+ assert(i1 <= i2);
+ #endif
+ // offset_ is 0-based, so no need for
+ // ( - TNT_BASE_OFFSET)
+ //
+ return const_Region1D<Array1D>(A_, i1+offset_,
+ offset_ + i2);
+ }
+
+
+ const_Region1D<Array1D> operator()(const Index1D &I)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(TNT_BASE_OFFSET<=I.lbound());
+ assert(I.ubound() <= dim() + (TNT_BASE_OFFSET-1));
+ assert(I.lbound() <= I.ubound());
+ #endif
+ return const_Region1D<Array1D>(A_, I.lbound()+offset_,
+ offset_ + I.ubound());
+ }
+
+
+ const T & operator() (Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(TNT_BASE_OFFSET <= i);
+ assert(i <= dim() + (TNT_BASE_OFFSET-1));
+ #endif
+ return A_(i+offset_);
+ }
+
+
+
+
+ };
+
+ template <class Array1D>
+ std::ostream& operator<<(std::ostream &s, const_Region1D<Array1D> &A)
+ {
+ Subscript N=A.dim();
+
+ for (Subscript i=1; i<=N; i++)
+ s << A(i) << endl;
+
+ return s;
+ }
+
+
+ } // namespace TNT
+
+ #endif
+ // const_Region1D_H
Index: tnt/region2d.h
===================================================================
RCS file: region2d.h
diff -N region2d.h
*** /dev/null Tue May 5 14:32:27 1998
--- region2d.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,467 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+ // 2D Regions for arrays and matrices
+
+ #ifndef REGION2D_H
+ #define REGION2D_H
+
+ #include "tnt/index.h"
+ #include <iostream>
+ #include <cassert>
+
+ namespace TNT
+ {
+
+ template <class Array2D>
+ class const_Region2D;
+
+
+ template <class Array2D>
+ class Region2D
+ {
+ protected:
+
+ Array2D & A_;
+ Subscript offset_[2]; // 0-offset internally
+ Subscript dim_[2];
+
+ public:
+ typedef typename Array2D::value_type T;
+ typedef Subscript size_type;
+ typedef T value_type;
+ typedef T element_type;
+ typedef T* pointer;
+ typedef T* iterator;
+ typedef T& reference;
+ typedef const T* const_iterator;
+ typedef const T& const_reference;
+
+ Array2D & array() { return A_; }
+ const Array2D & array() const { return A_; }
+ Subscript lbound() const { return A_.lbound(); }
+ Subscript num_rows() const { return dim_[0]; }
+ Subscript num_cols() const { return dim_[1]; }
+ Subscript offset(Subscript i) const // 1-offset
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( A_.lbound() <= i);
+ assert( i<= dim_[0] + A_.lbound()-1);
+ #endif
+ return offset_[i-A_.lbound()];
+ }
+
+ Subscript dim(Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( A_.lbound() <= i);
+ assert( i<= dim_[0] + A_.lbound()-1);
+ #endif
+ return dim_[i-A_.lbound()];
+ }
+
+
+
+ Region2D(Array2D &A, Subscript i1, Subscript i2, Subscript j1,
+ Subscript j2) : A_(A)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( i1 <= i2 );
+ assert( j1 <= j2);
+ assert( A.lbound() <= i1);
+ assert( i2<= A.dim(A.lbound()) + A.lbound()-1);
+ assert( A.lbound() <= j1);
+ assert( j2<= A.dim(A.lbound()+1) + A.lbound()-1 );
+ #endif
+
+
+ offset_[0] = i1-A.lbound();
+ offset_[1] = j1-A.lbound();
+ dim_[0] = i2-i1+1;
+ dim_[1] = j2-j1+1;
+ }
+
+ Region2D(Array2D &A, const Index1D &I, const Index1D &J) : A_(A)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( I.lbound() <= I.ubound() );
+ assert( J.lbound() <= J.ubound() );
+ assert( A.lbound() <= I.lbound());
+ assert( I.ubound()<= A.dim(A.lbound()) + A.lbound()-1);
+ assert( A.lbound() <= J.lbound());
+ assert( J.ubound() <= A.dim(A.lbound()+1) + A.lbound()-1 );
+ #endif
+
+ offset_[0] = I.lbound()-A.lbound();
+ offset_[1] = J.lbound()-A.lbound();
+ dim_[0] = I.ubound() - I.lbound() + 1;
+ dim_[1] = J.ubound() - J.lbound() + 1;
+ }
+
+ Region2D(Region2D<Array2D> &A, Subscript i1, Subscript i2,
+ Subscript j1, Subscript j2) : A_(A.A_)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( i1 <= i2 );
+ assert( j1 <= j2);
+ assert( A.lbound() <= i1);
+ assert( i2<= A.dim(A.lbound()) + A.lbound()-1);
+ assert( A.lbound() <= j1);
+ assert( j2<= A.dim(A.lbound()+1) + A.lbound()-1 );
+ #endif
+ offset_[0] = (i1 - A.lbound()) + A.offset_[0];
+ offset_[1] = (j1 - A.lbound()) + A.offset_[1];
+ dim_[0] = i2-i1 + 1;
+ dim_[1] = j2-j1+1;
+ }
+
+ Region2D<Array2D> operator()(Subscript i1, Subscript i2,
+ Subscript j1, Subscript j2)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( i1 <= i2 );
+ assert( j1 <= j2);
+ assert( A_.lbound() <= i1);
+ assert( i2<= dim_[0] + A_.lbound()-1);
+ assert( A_.lbound() <= j1);
+ assert( j2<= dim_[1] + A_.lbound()-1 );
+ #endif
+ return Region2D<Array2D>(A_,
+ i1+offset_[0], offset_[0] + i2,
+ j1+offset_[1], offset_[1] + j2);
+ }
+
+
+ Region2D<Array2D> operator()(const Index1D &I,
+ const Index1D &J)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( I.lbound() <= I.ubound() );
+ assert( J.lbound() <= J.ubound() );
+ assert( A_.lbound() <= I.lbound());
+ assert( I.ubound()<= dim_[0] + A_.lbound()-1);
+ assert( A_.lbound() <= J.lbound());
+ assert( J.ubound() <= dim_[1] + A_.lbound()-1 );
+ #endif
+
+ return Region2D<Array2D>(A_, I.lbound()+offset_[0],
+ offset_[0] + I.ubound(), offset_[1]+J.lbound(),
+ offset_[1] + J.ubound());
+ }
+
+ inline T & operator()(Subscript i, Subscript j)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( A_.lbound() <= i);
+ assert( i<= dim_[0] + A_.lbound()-1);
+ assert( A_.lbound() <= j);
+ assert( j<= dim_[1] + A_.lbound()-1 );
+ #endif
+ return A_(i+offset_[0], j+offset_[1]);
+ }
+
+ inline const T & operator() (Subscript i, Subscript j) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( A_.lbound() <= i);
+ assert( i<= dim_[0] + A_.lbound()-1);
+ assert( A_.lbound() <= j);
+ assert( j<= dim_[1] + A_.lbound()-1 );
+ #endif
+ return A_(i+offset_[0], j+offset_[1]);
+ }
+
+
+ Region2D<Array2D> & operator=(const Region2D<Array2D> &R)
+ {
+ Subscript M = num_rows();
+ Subscript N = num_cols();
+
+ // make sure both sides conform
+ assert(M == R.num_rows());
+ assert(N == R.num_cols());
+
+
+ Subscript start = R.lbound();
+ Subscript Mend = start + M - 1;
+ Subscript Nend = start + N - 1;
+ for (Subscript i=start; i<=Mend; i++)
+ for (Subscript j=start; j<=Nend; j++)
+ (*this)(i,j) = R(i,j);
+
+ return *this;
+ }
+
+ Region2D<Array2D> & operator=(const const_Region2D<Array2D> &R)
+ {
+ Subscript M = num_rows();
+ Subscript N = num_cols();
+
+ // make sure both sides conform
+ assert(M == R.num_rows());
+ assert(N == R.num_cols());
+
+
+ Subscript start = R.lbound();
+ Subscript Mend = start + M - 1;
+ Subscript Nend = start + N - 1;
+ for (Subscript i=start; i<=Mend; i++)
+ for (Subscript j=start; j<=Nend; j++)
+ (*this)(i,j) = R(i,j);
+
+ return *this;
+ }
+
+ Region2D<Array2D> & operator=(const Array2D &R)
+ {
+ Subscript M = num_rows();
+ Subscript N = num_cols();
+
+ // make sure both sides conform
+ assert(M == R.num_rows());
+ assert(N == R.num_cols());
+
+
+ Subscript start = R.lbound();
+ Subscript Mend = start + M - 1;
+ Subscript Nend = start + N - 1;
+ for (Subscript i=start; i<=Mend; i++)
+ for (Subscript j=start; j<=Nend; j++)
+ (*this)(i,j) = R(i,j);
+
+ return *this;
+ }
+
+ Region2D<Array2D> & operator=(const T &scalar)
+ {
+ Subscript start = lbound();
+ Subscript Mend = lbound() + num_rows() - 1;
+ Subscript Nend = lbound() + num_cols() - 1;
+
+
+ for (Subscript i=start; i<=Mend; i++)
+ for (Subscript j=start; j<=Nend; j++)
+ (*this)(i,j) = scalar;
+
+ return *this;
+ }
+
+ };
+
+ //************************
+
+ template <class Array2D>
+ class const_Region2D
+ {
+ protected:
+
+ const Array2D & A_;
+ Subscript offset_[2]; // 0-offset internally
+ Subscript dim_[2];
+
+ public:
+ typedef typename Array2D::value_type T;
+ typedef T value_type;
+ typedef T element_type;
+ typedef const T* const_iterator;
+ typedef const T& const_reference;
+
+ const Array2D & array() const { return A_; }
+ Subscript lbound() const { return A_.lbound(); }
+ Subscript num_rows() const { return dim_[0]; }
+ Subscript num_cols() const { return dim_[1]; }
+ Subscript offset(Subscript i) const // 1-offset
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( TNT_BASE_OFFSET <= i);
+ assert( i<= dim_[0] + TNT_BASE_OFFSET-1);
+ #endif
+ return offset_[i-TNT_BASE_OFFSET];
+ }
+
+ Subscript dim(Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( TNT_BASE_OFFSET <= i);
+ assert( i<= dim_[0] + TNT_BASE_OFFSET-1);
+ #endif
+ return dim_[i-TNT_BASE_OFFSET];
+ }
+
+
+ const_Region2D(const Array2D &A, Subscript i1, Subscript i2,
+ Subscript j1, Subscript j2) : A_(A)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( i1 <= i2 );
+ assert( j1 <= j2);
+ assert( TNT_BASE_OFFSET <= i1);
+ assert( i2<= A.dim(TNT_BASE_OFFSET) + TNT_BASE_OFFSET-1);
+ assert( TNT_BASE_OFFSET <= j1);
+ assert( j2<= A.dim(TNT_BASE_OFFSET+1) + TNT_BASE_OFFSET-1 );
+ #endif
+
+ offset_[0] = i1-TNT_BASE_OFFSET;
+ offset_[1] = j1-TNT_BASE_OFFSET;
+ dim_[0] = i2-i1+1;
+ dim_[1] = j2-j1+1;
+ }
+
+ const_Region2D(const Array2D &A, const Index1D &I, const Index1D &J)
+ : A_(A)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( I.lbound() <= I.ubound() );
+ assert( J.lbound() <= J.ubound() );
+ assert( TNT_BASE_OFFSET <= I.lbound());
+ assert( I.ubound()<= A.dim(TNT_BASE_OFFSET) + TNT_BASE_OFFSET-1);
+ assert( TNT_BASE_OFFSET <= J.lbound());
+ assert( J.ubound() <= A.dim(TNT_BASE_OFFSET+1) + TNT_BASE_OFFSET-1 );
+ #endif
+
+ offset_[0] = I.lbound()-TNT_BASE_OFFSET;
+ offset_[1] = J.lbound()-TNT_BASE_OFFSET;
+ dim_[0] = I.ubound() - I.lbound() + 1;
+ dim_[1] = J.ubound() - J.lbound() + 1;
+ }
+
+
+ const_Region2D(const_Region2D<Array2D> &A, Subscript i1,
+ Subscript i2,
+ Subscript j1, Subscript j2) : A_(A.A_)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( i1 <= i2 );
+ assert( j1 <= j2);
+ assert( TNT_BASE_OFFSET <= i1);
+ assert( i2<= A.dim(TNT_BASE_OFFSET) + TNT_BASE_OFFSET-1);
+ assert( TNT_BASE_OFFSET <= j1);
+ assert( j2<= A.dim(TNT_BASE_OFFSET+1) + TNT_BASE_OFFSET-1 );
+ #endif
+ offset_[0] = (i1 - TNT_BASE_OFFSET) + A.offset_[0];
+ offset_[1] = (j1 - TNT_BASE_OFFSET) + A.offset_[1];
+ dim_[0] = i2-i1 + 1;
+ dim_[1] = j2-j1+1;
+ }
+
+ const_Region2D<Array2D> operator()(Subscript i1, Subscript i2,
+ Subscript j1, Subscript j2)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( i1 <= i2 );
+ assert( j1 <= j2);
+ assert( TNT_BASE_OFFSET <= i1);
+ assert( i2<= dim_[0] + TNT_BASE_OFFSET-1);
+ assert( TNT_BASE_OFFSET <= j1);
+ assert( j2<= dim_[0] + TNT_BASE_OFFSET-1 );
+ #endif
+ return const_Region2D<Array2D>(A_,
+ i1+offset_[0], offset_[0] + i2,
+ j1+offset_[1], offset_[1] + j2);
+ }
+
+
+ const_Region2D<Array2D> operator()(const Index1D &I,
+ const Index1D &J)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( I.lbound() <= I.ubound() );
+ assert( J.lbound() <= J.ubound() );
+ assert( TNT_BASE_OFFSET <= I.lbound());
+ assert( I.ubound()<= dim_[0] + TNT_BASE_OFFSET-1);
+ assert( TNT_BASE_OFFSET <= J.lbound());
+ assert( J.ubound() <= dim_[1] + TNT_BASE_OFFSET-1 );
+ #endif
+
+ return const_Region2D<Array2D>(A_, I.lbound()+offset_[0],
+ offset_[0] + I.ubound(), offset_[1]+J.lbound(),
+ offset_[1] + J.ubound());
+ }
+
+
+ inline const T & operator() (Subscript i, Subscript j) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( TNT_BASE_OFFSET <= i);
+ assert( i<= dim_[0] + TNT_BASE_OFFSET-1);
+ assert( TNT_BASE_OFFSET <= j);
+ assert( j<= dim_[1] + TNT_BASE_OFFSET-1 );
+ #endif
+ return A_(i+offset_[0], j+offset_[1]);
+ }
+
+ };
+
+
+ // ************** std::ostream algorithms *******************************
+
+ template <class Array2D>
+ std::ostream& operator<<(std::ostream &s, const const_Region2D<Array2D> &A)
+ {
+ Subscript start = A.lbound();
+ Subscript Mend=A.lbound()+ A.num_rows() - 1;
+ Subscript Nend=A.lbound() + A.num_cols() - 1;
+
+
+ s << A.num_rows() << " " << A.num_cols() << "\n";
+ for (Subscript i=start; i<=Mend; i++)
+ {
+ for (Subscript j=start; j<=Nend; j++)
+ {
+ s << A(i,j) << " ";
+ }
+ s << "\n";
+ }
+
+
+ return s;
+ }
+
+ template <class Array2D>
+ std::ostream& operator<<(std::ostream &s, const Region2D<Array2D> &A)
+ {
+ Subscript start = A.lbound();
+ Subscript Mend=A.lbound()+ A.num_rows() - 1;
+ Subscript Nend=A.lbound() + A.num_cols() - 1;
+
+
+ s << A.num_rows() << " " << A.num_cols() << "\n";
+ for (Subscript i=start; i<=Mend; i++)
+ {
+ for (Subscript j=start; j<=Nend; j++)
+ {
+ s << A(i,j) << " ";
+ }
+ s << "\n";
+ }
+
+
+ return s;
+
+ }
+
+ } // namespace TNT
+
+ #endif
+ // REGION2D_H
Index: tnt/stopwatch.h
===================================================================
RCS file: stopwatch.h
diff -N stopwatch.h
*** /dev/null Tue May 5 14:32:27 1998
--- stopwatch.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,83 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ #ifndef STPWATCH_H
+ #define STPWATCH_H
+
+ // for clock() and CLOCKS_PER_SEC
+ #include <ctime>
+
+ namespace TNT
+ {
+
+ /* Simple stopwatch object:
+
+ void start() : start timing
+ double stop() : stop timing
+ void reset() : set elapsed time to 0.0
+ double read() : read elapsed time (in seconds)
+
+ */
+
+ inline double seconds(void)
+ {
+ static const double secs_per_tick = 1.0 / CLOCKS_PER_SEC;
+ return ( (double) clock() ) * secs_per_tick;
+ }
+
+
+ class stopwatch {
+ private:
+ int running;
+ double last_time;
+ double total;
+
+ public:
+ stopwatch() : running(0), last_time(0.0), total(0.0) {}
+ void reset() { running = 0; last_time = 0.0; total=0.0; }
+ void start() { if (!running) { last_time = seconds(); running = 1;}}
+ double stop() { if (running)
+ {
+ total += seconds() - last_time;
+ running = 0;
+ }
+ return total;
+ }
+ double read() { if (running)
+ {
+ total+= seconds() - last_time;
+ last_time = seconds();
+ }
+ return total;
+ }
+
+ };
+
+ } // namespace TNT
+
+ #endif
+
+
+
Index: tnt/subscript.h
===================================================================
RCS file: subscript.h
diff -N subscript.h
*** /dev/null Tue May 5 14:32:27 1998
--- subscript.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,58 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+ #ifndef SUBSCRPT_H
+ #define SUBSCRPT_H
+
+
+ //---------------------------------------------------------------------
+ // This definition describes the default TNT data type used for
+ // indexing into TNT matrices and vectors. The data type should
+ // be wide enough to index into large arrays. It defaults to an
+ // "int", but can be overriden at compile time redefining TNT_SUBSCRIPT_TYPE,
+ // e.g.
+ //
+ // g++ -DTNT_SUBSCRIPT_TYPE='unsigned int' ...
+ //
+ //---------------------------------------------------------------------
+ //
+
+ #ifndef TNT_SUBSCRIPT_TYPE
+ #define TNT_SUBSCRIPT_TYPE int
+ #endif
+
+ namespace TNT
+ {
+ typedef TNT_SUBSCRIPT_TYPE Subscript;
+ }
+
+
+ // () indexing in TNT means 1-offset, i.e. x(1) and A(1,1) are the
+ // first elements. This offset is left as a macro for future
+ // purposes, but should not be changed in the current release.
+ //
+ //
+ #define TNT_BASE_OFFSET (1)
+
+ #endif
Index: tnt/tnt.h
===================================================================
RCS file: tnt.h
diff -N tnt.h
*** /dev/null Tue May 5 14:32:27 1998
--- tnt.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,90 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+ #ifndef TNT_H
+ #define TNT_H
+
+ //---------------------------------------------------------------------
+ // tnt.h TNT general header file. Defines default types
+ // and conventions.
+ //---------------------------------------------------------------------
+
+ //---------------------------------------------------------------------
+ // Include current version
+ //---------------------------------------------------------------------
+ #include "tnt/version.h"
+
+ //---------------------------------------------------------------------
+ // Define the data type used for matrix and vector Subscripts.
+ // This will default to "int", but it can be overridden at compile time,
+ // e.g.
+ //
+ // g++ -DTNT_SUBSCRIPT_TYPE='unsigned long' ...
+ //
+ // See subscript.h for details.
+ //---------------------------------------------------------------------
+
+ #include "tnt/subscript.h"
+
+
+
+ //---------------------------------------------------------------------
+ // Define this macro if you want TNT to ensure all refernces
+ // are within the bounds of the array. This encurs a run-time
+ // overhead, of course, but is recommended while developing
+ // code. It can be turned off for production runs.
+ //
+ // #define TNT_BOUNDS_CHECK
+ //---------------------------------------------------------------------
+ //
+ #define TNT_BOUNDS_CHECK
+ #ifdef TNT_NO_BOUNDS_CHECK
+ #undef TNT_BOUNDS_CHECK
+ #endif
+
+ //---------------------------------------------------------------------
+ // Define this macro if you want to utilize matrix and vector
+ // regions. This is typically on, but you can save some
+ // compilation time by turning it off. If you do this and
+ // attempt to use regions you will get an error message.
+ //
+ // #define TNT_USE_REGIONS
+ //---------------------------------------------------------------------
+ //
+ #define TNT_USE_REGIONS
+
+ //---------------------------------------------------------------------
+ //
+ //---------------------------------------------------------------------
+ // If your system doesn't have abs(), min(), and max() uncomment the following.
+ //---------------------------------------------------------------------
+ //
+ //
+ //#define __NEED_ABS_MIN_MAX_
+
+ #include "tnt/tntmath.h"
+
+
+ #endif
+ // TNT_H
Index: tnt/tnt.h.patch
===================================================================
RCS file: tnt.h.patch
diff -N tnt.h.patch
*** /dev/null Tue May 5 14:32:27 1998
--- tnt.h.patch Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,47 ----
+ *** tnt.h.~1~ Sun Aug 6 21:52:04 2000
+ --- tnt.h Tue Aug 7 19:27:13 2001
+ ***************
+ *** 38,45 ****
+ //---------------------------------------------------------------------
+ // Define the data type used for matrix and vector Subscripts.
+ ! // This will default to "int", but it can be overriden at compile time,
+ // e.g.
+ //
+ ! // g++ -DTNT_SUBSCRIPT_TYPE='unsinged long' ...
+ //
+ // See subscript.h for details.
+ --- 38,45 ----
+ //---------------------------------------------------------------------
+ // Define the data type used for matrix and vector Subscripts.
+ ! // This will default to "int", but it can be overridden at compile time,
+ // e.g.
+ //
+ ! // g++ -DTNT_SUBSCRIPT_TYPE='unsigned long' ...
+ //
+ // See subscript.h for details.
+ ***************
+ *** 51,55 ****
+
+ //---------------------------------------------------------------------
+ ! // Define this macro if you want TNT to ensure all refernces
+ // are within the bounds of the array. This encurs a run-time
+ // overhead, of course, but is recommended while developing
+ --- 51,55 ----
+
+ //---------------------------------------------------------------------
+ ! // Define this macro if you want TNT to ensure all refernces
+ // are within the bounds of the array. This encurs a run-time
+ // overhead, of course, but is recommended while developing
+ ***************
+ *** 78,82 ****
+ //
+ //---------------------------------------------------------------------
+ ! // if your system doesn't have abs() min(), and max() uncoment the following
+ //---------------------------------------------------------------------
+ //
+ --- 78,82 ----
+ //
+ //---------------------------------------------------------------------
+ ! // If your system doesn't have abs(), min(), and max() uncomment the following.
+ //---------------------------------------------------------------------
+ //
Index: tnt/tntmath.h
===================================================================
RCS file: tntmath.h
diff -N tntmath.h
*** /dev/null Tue May 5 14:32:27 1998
--- tntmath.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,85 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ // Header file for scalar math functions
+
+ #ifndef TNTMATH_H
+ #define TNTMATH_H
+
+ // conventional functions required by several matrix algorithms
+
+
+
+ namespace TNT
+ {
+
+ inline double abs(double t)
+ {
+ return ( t > 0 ? t : -t);
+ }
+
+ inline double min(double a, double b)
+ {
+ return (a < b ? a : b);
+ }
+
+ inline double max(double a, double b)
+ {
+ return (a > b ? a : b);
+ }
+
+ inline float abs(float t)
+ {
+ return ( t > 0 ? t : -t);
+ }
+
+ inline float min(float a, float b)
+ {
+ return (a < b ? a : b);
+ }
+
+ inline float max(float a, float b)
+ {
+ return (a > b ? a : b);
+ }
+
+ inline double sign(double a)
+ {
+ return (a > 0 ? 1.0 : -1.0);
+ }
+
+
+
+ inline float sign(float a)
+ {
+ return (a > 0.0 ? 1.0f : -1.0f);
+ }
+
+ } /* namespace TNT */
+
+
+ #endif
+ /* TNTMATH_H */
+
Index: tnt/tntreqs.h
===================================================================
RCS file: tntreqs.h
diff -N tntreqs.h
*** /dev/null Tue May 5 14:32:27 1998
--- tntreqs.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,70 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ // The requirements for a bare-bones vector class:
+ //
+ //
+ // o) must have 0-based [] indexing for const and
+ // non-const objects (i.e. operator[] defined)
+ //
+ // o) must have size() method to denote the number of
+ // elements
+ // o) must clean up after itself when destructed
+ // (i.e. no memory leaks)
+ //
+ // -) must have begin() and end() methods (The begin()
+ // method is necessary, because relying on
+ // &v_[0] may not work on a empty vector (i.e. v_ is NULL.)
+ //
+ // o) must be templated
+ // o) must have X::value_type defined to be the types of elements
+ // o) must have X::X(const &x) copy constructor (by *value*)
+ // o) must have X::X(int N) constructor to N-length vector
+ // (NOTE: this constructor need *NOT* initalize elements)
+ //
+ // -) must have X::X(int N, T scalar) constructor to initalize
+ // elements to value of "scalar".
+ //
+ // ( removed, because valarray<> class uses (scalar, N) rather
+ // than (N, scalar) )
+ // -) must have X::X(int N, const T* scalars) constructor to copy from
+ // any C linear array
+ //
+ // ( removed, because of same reverse order of valarray<> )
+ //
+ // o) must have assignment A=B, by value
+ //
+ // NOTE: this class is *NOT* meant to be derived from,
+ // so its methods (particularly indexing) need not be
+ // declared virtual.
+ //
+ //
+ // Some things it *DOES NOT* need to do are
+ //
+ // o) bounds checking
+ // o) array referencing (e.g. reference counting)
+ // o) support () indexing
+ // o) I/O
+ //
Index: tnt/transv.h
===================================================================
RCS file: transv.h
diff -N transv.h
*** /dev/null Tue May 5 14:32:27 1998
--- transv.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,160 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ // Matrix Transpose Views
+
+ #ifndef TRANSV_H
+ #define TRANSV_H
+
+ #include <iostream>
+ #include <cassert>
+ #include "tnt/vec.h"
+
+ namespace TNT
+ {
+
+ template <class Array2D>
+ class Transpose_View
+ {
+ protected:
+
+ const Array2D & A_;
+
+ public:
+
+ typedef typename Array2D::element_type T;
+ typedef T value_type;
+ typedef T element_type;
+ typedef T* pointer;
+ typedef T* iterator;
+ typedef T& reference;
+ typedef const T* const_iterator;
+ typedef const T& const_reference;
+
+
+ const Array2D & array() const { return A_; }
+ Subscript num_rows() const { return A_.num_cols();}
+ Subscript num_cols() const { return A_.num_rows();}
+ Subscript lbound() const { return A_.lbound(); }
+ Subscript dim(Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert( A_.lbound() <= i);
+ assert( i<= A_.lbound()+1);
+ #endif
+ if (i== A_.lbound())
+ return num_rows();
+ else
+ return num_cols();
+ }
+
+
+ Transpose_View(const Transpose_View<Array2D> &A) : A_(A.A_) {};
+ Transpose_View(const Array2D &A) : A_(A) {};
+
+
+ inline const typename Array2D::element_type & operator()(
+ Subscript i, Subscript j) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(lbound()<=i);
+ assert(i<=A_.num_cols() + lbound() - 1);
+ assert(lbound()<=j);
+ assert(j<=A_.num_rows() + lbound() - 1);
+ #endif
+
+ return A_(j,i);
+ }
+
+
+ };
+
+ template <class Matrix>
+ Transpose_View<Matrix> Transpose_view(const Matrix &A)
+ {
+ return Transpose_View<Matrix>(A);
+ }
+
+ template <class Matrix, class T>
+ Vector<T> matmult(
+ const Transpose_View<Matrix> & A,
+ const Vector<T> &B)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ assert(B.dim() == N);
+
+ Vector<T> x(N);
+
+ Subscript i, j;
+ T tmp = 0;
+
+ for (i=1; i<=M; i++)
+ {
+ tmp = 0;
+ for (j=1; j<=N; j++)
+ tmp += A(i,j) * B(j);
+ x(i) = tmp;
+ }
+
+ return x;
+ }
+
+ template <class Matrix, class T>
+ inline Vector<T> operator*(const Transpose_View<Matrix> & A, const Vector<T> &B)
+ {
+ return matmult(A,B);
+ }
+
+
+ template <class Matrix>
+ std::ostream& operator<<(std::ostream &s, const Transpose_View<Matrix> &A)
+ {
+ Subscript M=A.num_rows();
+ Subscript N=A.num_cols();
+
+ Subscript start = A.lbound();
+ Subscript Mend = M + A.lbound() - 1;
+ Subscript Nend = N + A.lbound() - 1;
+
+ s << M << " " << N << endl;
+ for (Subscript i=start; i<=Mend; i++)
+ {
+ for (Subscript j=start; j<=Nend; j++)
+ {
+ s << A(i,j) << " ";
+ }
+ s << endl;
+ }
+
+
+ return s;
+ }
+
+ } // namespace TNT
+
+ #endif
+ // TRANSV_H
Index: tnt/triang.h
===================================================================
RCS file: triang.h
diff -N triang.h
*** /dev/null Tue May 5 14:32:27 1998
--- triang.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,638 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ // Triangular Matrices (Views and Adpators)
+
+ #ifndef TRIANG_H
+ #define TRIANG_H
+
+ // default to use lower-triangular portions of arrays
+ // for symmetric matrices.
+
+ namespace TNT
+ {
+
+ template <class MaTRiX>
+ class LowerTriangularView
+ {
+ protected:
+
+
+ const MaTRiX &A_;
+ const typename MaTRiX::element_type zero_;
+
+ public:
+
+
+ typedef typename MaTRiX::const_reference const_reference;
+ typedef const typename MaTRiX::element_type element_type;
+ typedef const typename MaTRiX::element_type value_type;
+ typedef element_type T;
+
+ Subscript dim(Subscript d) const { return A_.dim(d); }
+ Subscript lbound() const { return A_.lbound(); }
+ Subscript num_rows() const { return A_.num_rows(); }
+ Subscript num_cols() const { return A_.num_cols(); }
+
+
+ // constructors
+
+ LowerTriangularView(/*const*/ MaTRiX &A) : A_(A), zero_(0) {}
+
+
+ inline const_reference get(Subscript i, Subscript j) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(lbound()<=i);
+ assert(i<=A_.num_rows() + lbound() - 1);
+ assert(lbound()<=j);
+ assert(j<=A_.num_cols() + lbound() - 1);
+ #endif
+ if (i<j)
+ return zero_;
+ else
+ return A_(i,j);
+ }
+
+
+ inline const_reference operator() (Subscript i, Subscript j) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(lbound()<=i);
+ assert(i<=A_.num_rows() + lbound() - 1);
+ assert(lbound()<=j);
+ assert(j<=A_.num_cols() + lbound() - 1);
+ #endif
+ if (i<j)
+ return zero_;
+ else
+ return A_(i,j);
+ }
+
+ #ifdef TNT_USE_REGIONS
+
+ typedef const_Region2D< LowerTriangularView<MaTRiX> >
+ const_Region;
+
+ const_Region operator()(/*const*/ Index1D &I,
+ /*const*/ Index1D &J) const
+ {
+ return const_Region(*this, I, J);
+ }
+
+ const_Region operator()(Subscript i1, Subscript i2,
+ Subscript j1, Subscript j2) const
+ {
+ return const_Region(*this, i1, i2, j1, j2);
+ }
+
+
+
+ #endif
+ // TNT_USE_REGIONS
+
+ };
+
+
+ /* *********** Lower_triangular_view() algorithms ****************** */
+
+ template <class MaTRiX, class VecToR>
+ VecToR matmult(/*const*/ LowerTriangularView<MaTRiX> &A, VecToR &x)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ assert(N == x.dim());
+
+ Subscript i, j;
+ typename MaTRiX::element_type sum=0.0;
+ VecToR result(M);
+
+ Subscript start = A.lbound();
+ Subscript Mend = M + A.lbound() -1 ;
+
+ for (i=start; i<=Mend; i++)
+ {
+ sum = 0.0;
+ for (j=start; j<=i; j++)
+ sum = sum + A(i,j)*x(j);
+ result(i) = sum;
+ }
+
+ return result;
+ }
+
+ template <class MaTRiX, class VecToR>
+ inline VecToR operator*(/*const*/ LowerTriangularView<MaTRiX> &A, VecToR &x)
+ {
+ return matmult(A,x);
+ }
+
+ template <class MaTRiX>
+ class UnitLowerTriangularView
+ {
+ protected:
+
+ const MaTRiX &A_;
+ const typename MaTRiX::element_type zero;
+ const typename MaTRiX::element_type one;
+
+ public:
+
+ typedef typename MaTRiX::const_reference const_reference;
+ typedef typename MaTRiX::element_type element_type;
+ typedef typename MaTRiX::element_type value_type;
+ typedef element_type T;
+
+ Subscript lbound() const { return 1; }
+ Subscript dim(Subscript d) const { return A_.dim(d); }
+ Subscript num_rows() const { return A_.num_rows(); }
+ Subscript num_cols() const { return A_.num_cols(); }
+
+
+ // constructors
+
+ UnitLowerTriangularView(/*const*/ MaTRiX &A) : A_(A), zero(0), one(1) {}
+
+
+ inline const_reference get(Subscript i, Subscript j) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i<=A_.dim(1));
+ assert(1<=j);
+ assert(j<=A_.dim(2));
+ assert(0<=i && i<A_.dim(0) && 0<=j && j<A_.dim(1));
+ #endif
+ if (i>j)
+ return A_(i,j);
+ else if (i==j)
+ return one;
+ else
+ return zero;
+ }
+
+
+ inline const_reference operator() (Subscript i, Subscript j) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i<=A_.dim(1));
+ assert(1<=j);
+ assert(j<=A_.dim(2));
+ #endif
+ if (i>j)
+ return A_(i,j);
+ else if (i==j)
+ return one;
+ else
+ return zero;
+ }
+
+
+ #ifdef TNT_USE_REGIONS
+ // These are the "index-aware" features
+
+ typedef const_Region2D< UnitLowerTriangularView<MaTRiX> >
+ const_Region;
+
+ const_Region operator()(/*const*/ Index1D &I,
+ /*const*/ Index1D &J) const
+ {
+ return const_Region(*this, I, J);
+ }
+
+ const_Region operator()(Subscript i1, Subscript i2,
+ Subscript j1, Subscript j2) const
+ {
+ return const_Region(*this, i1, i2, j1, j2);
+ }
+ #endif
+ // TNT_USE_REGIONS
+ };
+
+ template <class MaTRiX>
+ LowerTriangularView<MaTRiX> Lower_triangular_view(
+ /*const*/ MaTRiX &A)
+ {
+ return LowerTriangularView<MaTRiX>(A);
+ }
+
+
+ template <class MaTRiX>
+ UnitLowerTriangularView<MaTRiX> Unit_lower_triangular_view(
+ /*const*/ MaTRiX &A)
+ {
+ return UnitLowerTriangularView<MaTRiX>(A);
+ }
+
+ template <class MaTRiX, class VecToR>
+ VecToR matmult(/*const*/ UnitLowerTriangularView<MaTRiX> &A, VecToR &x)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ assert(N == x.dim());
+
+ Subscript i, j;
+ typename MaTRiX::element_type sum=0.0;
+ VecToR result(M);
+
+ Subscript start = A.lbound();
+ Subscript Mend = M + A.lbound() -1 ;
+
+ for (i=start; i<=Mend; i++)
+ {
+ sum = 0.0;
+ for (j=start; j<i; j++)
+ sum = sum + A(i,j)*x(j);
+ result(i) = sum + x(i);
+ }
+
+ return result;
+ }
+
+ template <class MaTRiX, class VecToR>
+ inline VecToR operator*(/*const*/ UnitLowerTriangularView<MaTRiX> &A, VecToR &x)
+ {
+ return matmult(A,x);
+ }
+
+
+ //********************** Algorithms *************************************
+
+
+
+ template <class MaTRiX>
+ std::ostream& operator<<(std::ostream &s, const LowerTriangularView<MaTRiX>&A)
+ {
+ Subscript M=A.num_rows();
+ Subscript N=A.num_cols();
+
+ s << M << " " << N << endl;
+
+ for (Subscript i=1; i<=M; i++)
+ {
+ for (Subscript j=1; j<=N; j++)
+ {
+ s << A(i,j) << " ";
+ }
+ s << endl;
+ }
+
+
+ return s;
+ }
+
+ template <class MaTRiX>
+ std::ostream& operator<<(std::ostream &s,
+ const UnitLowerTriangularView<MaTRiX>&A)
+ {
+ Subscript M=A.num_rows();
+ Subscript N=A.num_cols();
+
+ s << M << " " << N << endl;
+
+ for (Subscript i=1; i<=M; i++)
+ {
+ for (Subscript j=1; j<=N; j++)
+ {
+ s << A(i,j) << " ";
+ }
+ s << endl;
+ }
+
+
+ return s;
+ }
+
+
+
+ // ******************* Upper Triangular Section **************************
+
+ template <class MaTRiX>
+ class UpperTriangularView
+ {
+ protected:
+
+
+ /*const*/ MaTRiX &A_;
+ /*const*/ typename MaTRiX::element_type zero_;
+
+ public:
+
+
+ typedef typename MaTRiX::const_reference const_reference;
+ typedef /*const*/ typename MaTRiX::element_type element_type;
+ typedef /*const*/ typename MaTRiX::element_type value_type;
+ typedef element_type T;
+
+ Subscript dim(Subscript d) const { return A_.dim(d); }
+ Subscript lbound() const { return A_.lbound(); }
+ Subscript num_rows() const { return A_.num_rows(); }
+ Subscript num_cols() const { return A_.num_cols(); }
+
+
+ // constructors
+
+ UpperTriangularView(/*const*/ MaTRiX &A) : A_(A), zero_(0) {}
+
+
+ inline const_reference get(Subscript i, Subscript j) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(lbound()<=i);
+ assert(i<=A_.num_rows() + lbound() - 1);
+ assert(lbound()<=j);
+ assert(j<=A_.num_cols() + lbound() - 1);
+ #endif
+ if (i>j)
+ return zero_;
+ else
+ return A_(i,j);
+ }
+
+
+ inline const_reference operator() (Subscript i, Subscript j) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(lbound()<=i);
+ assert(i<=A_.num_rows() + lbound() - 1);
+ assert(lbound()<=j);
+ assert(j<=A_.num_cols() + lbound() - 1);
+ #endif
+ if (i>j)
+ return zero_;
+ else
+ return A_(i,j);
+ }
+
+ #ifdef TNT_USE_REGIONS
+
+ typedef const_Region2D< UpperTriangularView<MaTRiX> >
+ const_Region;
+
+ const_Region operator()(const Index1D &I,
+ const Index1D &J) const
+ {
+ return const_Region(*this, I, J);
+ }
+
+ const_Region operator()(Subscript i1, Subscript i2,
+ Subscript j1, Subscript j2) const
+ {
+ return const_Region(*this, i1, i2, j1, j2);
+ }
+
+
+
+ #endif
+ // TNT_USE_REGIONS
+
+ };
+
+
+ /* *********** Upper_triangular_view() algorithms ****************** */
+
+ template <class MaTRiX, class VecToR>
+ VecToR matmult(/*const*/ UpperTriangularView<MaTRiX> &A, VecToR &x)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ assert(N == x.dim());
+
+ Subscript i, j;
+ typename VecToR::element_type sum=0.0;
+ VecToR result(M);
+
+ Subscript start = A.lbound();
+ Subscript Mend = M + A.lbound() -1 ;
+
+ for (i=start; i<=Mend; i++)
+ {
+ sum = 0.0;
+ for (j=i; j<=N; j++)
+ sum = sum + A(i,j)*x(j);
+ result(i) = sum;
+ }
+
+ return result;
+ }
+
+ template <class MaTRiX, class VecToR>
+ inline VecToR operator*(/*const*/ UpperTriangularView<MaTRiX> &A, VecToR &x)
+ {
+ return matmult(A,x);
+ }
+
+ template <class MaTRiX>
+ class UnitUpperTriangularView
+ {
+ protected:
+
+ const MaTRiX &A_;
+ const typename MaTRiX::element_type zero;
+ const typename MaTRiX::element_type one;
+
+ public:
+
+ typedef typename MaTRiX::const_reference const_reference;
+ typedef typename MaTRiX::element_type element_type;
+ typedef typename MaTRiX::element_type value_type;
+ typedef element_type T;
+
+ Subscript lbound() const { return 1; }
+ Subscript dim(Subscript d) const { return A_.dim(d); }
+ Subscript num_rows() const { return A_.num_rows(); }
+ Subscript num_cols() const { return A_.num_cols(); }
+
+
+ // constructors
+
+ UnitUpperTriangularView(/*const*/ MaTRiX &A) : A_(A), zero(0), one(1) {}
+
+
+ inline const_reference get(Subscript i, Subscript j) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i<=A_.dim(1));
+ assert(1<=j);
+ assert(j<=A_.dim(2));
+ assert(0<=i && i<A_.dim(0) && 0<=j && j<A_.dim(1));
+ #endif
+ if (i<j)
+ return A_(i,j);
+ else if (i==j)
+ return one;
+ else
+ return zero;
+ }
+
+
+ inline const_reference operator() (Subscript i, Subscript j) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i<=A_.dim(1));
+ assert(1<=j);
+ assert(j<=A_.dim(2));
+ #endif
+ if (i<j)
+ return A_(i,j);
+ else if (i==j)
+ return one;
+ else
+ return zero;
+ }
+
+
+ #ifdef TNT_USE_REGIONS
+ // These are the "index-aware" features
+
+ typedef const_Region2D< UnitUpperTriangularView<MaTRiX> >
+ const_Region;
+
+ const_Region operator()(const Index1D &I,
+ const Index1D &J) const
+ {
+ return const_Region(*this, I, J);
+ }
+
+ const_Region operator()(Subscript i1, Subscript i2,
+ Subscript j1, Subscript j2) const
+ {
+ return const_Region(*this, i1, i2, j1, j2);
+ }
+ #endif
+ // TNT_USE_REGIONS
+ };
+
+ template <class MaTRiX>
+ UpperTriangularView<MaTRiX> Upper_triangular_view(
+ /*const*/ MaTRiX &A)
+ {
+ return UpperTriangularView<MaTRiX>(A);
+ }
+
+
+ template <class MaTRiX>
+ UnitUpperTriangularView<MaTRiX> Unit_upper_triangular_view(
+ /*const*/ MaTRiX &A)
+ {
+ return UnitUpperTriangularView<MaTRiX>(A);
+ }
+
+ template <class MaTRiX, class VecToR>
+ VecToR matmult(/*const*/ UnitUpperTriangularView<MaTRiX> &A, VecToR &x)
+ {
+ Subscript M = A.num_rows();
+ Subscript N = A.num_cols();
+
+ assert(N == x.dim());
+
+ Subscript i, j;
+ typename VecToR::element_type sum=0.0;
+ VecToR result(M);
+
+ Subscript start = A.lbound();
+ Subscript Mend = M + A.lbound() -1 ;
+
+ for (i=start; i<=Mend; i++)
+ {
+ sum = x(i);
+ for (j=i+1; j<=N; j++)
+ sum = sum + A(i,j)*x(j);
+ result(i) = sum + x(i);
+ }
+
+ return result;
+ }
+
+ template <class MaTRiX, class VecToR>
+ inline VecToR operator*(/*const*/ UnitUpperTriangularView<MaTRiX> &A, VecToR &x)
+ {
+ return matmult(A,x);
+ }
+
+
+ //********************** Algorithms *************************************
+
+
+
+ template <class MaTRiX>
+ std::ostream& operator<<(std::ostream &s,
+ /*const*/ UpperTriangularView<MaTRiX>&A)
+ {
+ Subscript M=A.num_rows();
+ Subscript N=A.num_cols();
+
+ s << M << " " << N << endl;
+
+ for (Subscript i=1; i<=M; i++)
+ {
+ for (Subscript j=1; j<=N; j++)
+ {
+ s << A(i,j) << " ";
+ }
+ s << endl;
+ }
+
+
+ return s;
+ }
+
+ template <class MaTRiX>
+ std::ostream& operator<<(std::ostream &s,
+ /*const*/ UnitUpperTriangularView<MaTRiX>&A)
+ {
+ Subscript M=A.num_rows();
+ Subscript N=A.num_cols();
+
+ s << M << " " << N << endl;
+
+ for (Subscript i=1; i<=M; i++)
+ {
+ for (Subscript j=1; j<=N; j++)
+ {
+ s << A(i,j) << " ";
+ }
+ s << endl;
+ }
+
+
+ return s;
+ }
+
+ } // namespace TNT
+
+
+
+
+
+ #endif
+ //TRIANG_H
+
Index: tnt/trisolve.h
===================================================================
RCS file: trisolve.h
diff -N trisolve.h
*** /dev/null Tue May 5 14:32:27 1998
--- trisolve.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,184 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ // Triangular Solves
+
+ #ifndef TRISLV_H
+ #define TRISLV_H
+
+
+ #include "triang.h"
+
+ namespace TNT
+ {
+
+ template <class MaTriX, class VecToR>
+ VecToR Lower_triangular_solve(/*const*/ MaTriX &A, /*const*/ VecToR &b)
+ {
+ Subscript N = A.num_rows();
+
+ // make sure matrix sizes agree; A must be square
+
+ assert(A.num_cols() == N);
+ assert(b.dim() == N);
+
+ VecToR x(N);
+
+ Subscript i;
+ for (i=1; i<=N; i++)
+ {
+ typename MaTriX::element_type tmp=0;
+
+ for (Subscript j=1; j<i; j++)
+ tmp = tmp + A(i,j)*x(j);
+
+ x(i) = (b(i) - tmp)/ A(i,i);
+ }
+
+ return x;
+ }
+
+
+ template <class MaTriX, class VecToR>
+ VecToR Unit_lower_triangular_solve(/*const*/ MaTriX &A, /*const*/ VecToR &b)
+ {
+ Subscript N = A.num_rows();
+
+ // make sure matrix sizes agree; A must be square
+
+ assert(A.num_cols() == N);
+ assert(b.dim() == N);
+
+ VecToR x(N);
+
+ Subscript i;
+ for (i=1; i<=N; i++)
+ {
+
+ typename MaTriX::element_type tmp=0;
+
+ for (Subscript j=1; j<i; j++)
+ tmp = tmp + A(i,j)*x(j);
+
+ x(i) = b(i) - tmp;
+ }
+
+ return x;
+ }
+
+
+ template <class MaTriX, class VecToR>
+ VecToR linear_solve(/*const*/ LowerTriangularView<MaTriX> &A,
+ /*const*/ VecToR &b)
+ {
+ return Lower_triangular_solve(A, b);
+ }
+
+ template <class MaTriX, class VecToR>
+ VecToR linear_solve(/*const*/ UnitLowerTriangularView<MaTriX> &A,
+ /*const*/ VecToR &b)
+ {
+ return Unit_lower_triangular_solve(A, b);
+ }
+
+
+
+ //********************** Upper triangular section ****************
+
+ template <class MaTriX, class VecToR>
+ VecToR Upper_triangular_solve(/*const*/ MaTriX &A, /*const*/ VecToR &b)
+ {
+ Subscript N = A.num_rows();
+
+ // make sure matrix sizes agree; A must be square
+
+ assert(A.num_cols() == N);
+ assert(b.dim() == N);
+
+ VecToR x(N);
+
+ Subscript i;
+ for (i=N; i>=1; i--)
+ {
+
+ typename MaTriX::element_type tmp=0;
+
+ for (Subscript j=i+1; j<=N; j++)
+ tmp = tmp + A(i,j)*x(j);
+
+ x(i) = (b(i) - tmp)/ A(i,i);
+ }
+
+ return x;
+ }
+
+
+ template <class MaTriX, class VecToR>
+ VecToR Unit_upper_triangular_solve(/*const*/ MaTriX &A, /*const*/ VecToR &b)
+ {
+ Subscript N = A.num_rows();
+
+ // make sure matrix sizes agree; A must be square
+
+ assert(A.num_cols() == N);
+ assert(b.dim() == N);
+
+ VecToR x(N);
+
+ Subscript i;
+ for (i=N; i>=1; i--)
+ {
+
+ typename MaTriX::element_type tmp=0;
+
+ for (Subscript j=i+1; j<i; j++)
+ tmp = tmp + A(i,j)*x(j);
+
+ x(i) = b(i) - tmp;
+ }
+
+ return x;
+ }
+
+
+ template <class MaTriX, class VecToR>
+ VecToR linear_solve(/*const*/ UpperTriangularView<MaTriX> &A,
+ /*const*/ VecToR &b)
+ {
+ return Upper_triangular_solve(A, b);
+ }
+
+ template <class MaTriX, class VecToR>
+ VecToR linear_solve(/*const*/ UnitUpperTriangularView<MaTriX> &A,
+ /*const*/ VecToR &b)
+ {
+ return Unit_upper_triangular_solve(A, b);
+ }
+
+
+ } // namespace TNT
+
+ #endif
+ // TRISLV_H
Index: tnt/vec.h
===================================================================
RCS file: vec.h
diff -N vec.h
*** /dev/null Tue May 5 14:32:27 1998
--- vec.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,403 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+ // Basic TNT numerical vector (0-based [i] AND 1-based (i) indexing )
+ //
+
+ #ifndef VEC_H
+ #define VEC_H
+
+ #include "tnt/subscript.h"
+ #include <cstdlib>
+ #include <cassert>
+ #include <iostream>
+ #include <sstream>
+
+ namespace TNT
+ {
+
+ template <class T>
+ class Vector
+ {
+
+
+ public:
+
+ typedef Subscript size_type;
+ typedef T value_type;
+ typedef T element_type;
+ typedef T* pointer;
+ typedef T* iterator;
+ typedef T& reference;
+ typedef const T* const_iterator;
+ typedef const T& const_reference;
+
+ Subscript lbound() const { return 1;}
+
+ protected:
+ T* v_;
+ T* vm1_; // pointer adjustment for optimized 1-offset indexing
+ Subscript n_;
+
+ // internal helper function to create the array
+ // of row pointers
+
+ void initialize(Subscript N)
+ {
+ // adjust pointers so that they are 1-offset:
+ // v_[] is the internal contiguous array, it is still 0-offset
+ //
+ assert(v_ == NULL);
+ v_ = new T[N];
+ assert(v_ != NULL);
+ vm1_ = v_-1;
+ n_ = N;
+ }
+
+ void copy(const T* v)
+ {
+ Subscript N = n_;
+ Subscript i;
+
+ #ifdef TNT_UNROLL_LOOPS
+ Subscript Nmod4 = N & 3;
+ Subscript N4 = N - Nmod4;
+
+ for (i=0; i<N4; i+=4)
+ {
+ v_[i] = v[i];
+ v_[i+1] = v[i+1];
+ v_[i+2] = v[i+2];
+ v_[i+3] = v[i+3];
+ }
+
+ for (i=N4; i< N; i++)
+ v_[i] = v[i];
+ #else
+
+ for (i=0; i< N; i++)
+ v_[i] = v[i];
+ #endif
+ }
+
+ void set(const T& val)
+ {
+ Subscript N = n_;
+ Subscript i;
+
+ #ifdef TNT_UNROLL_LOOPS
+ Subscript Nmod4 = N & 3;
+ Subscript N4 = N - Nmod4;
+
+ for (i=0; i<N4; i+=4)
+ {
+ v_[i] = val;
+ v_[i+1] = val;
+ v_[i+2] = val;
+ v_[i+3] = val;
+ }
+
+ for (i=N4; i< N; i++)
+ v_[i] = val;
+ #else
+
+ for (i=0; i< N; i++)
+ v_[i] = val;
+
+ #endif
+ }
+
+
+
+ void destroy()
+ {
+ /* do nothing, if no memory has been previously allocated */
+ if (v_ == NULL) return ;
+
+ /* if we are here, then matrix was previously allocated */
+ delete [] (v_);
+
+ v_ = NULL;
+ vm1_ = NULL;
+ }
+
+
+ public:
+
+ // access
+
+ iterator begin() { return v_;}
+ iterator end() { return v_ + n_; }
+ const iterator begin() const { return v_;}
+ const iterator end() const { return v_ + n_; }
+
+ // destructor
+
+ ~Vector()
+ {
+ destroy();
+ }
+
+ // constructors
+
+ Vector() : v_(0), vm1_(0), n_(0) {};
+
+ Vector(const Vector<T> &A) : v_(0), vm1_(0), n_(0)
+ {
+ initialize(A.n_);
+ copy(A.v_);
+ }
+
+ Vector(Subscript N, const T& value = T()) : v_(0), vm1_(0), n_(0)
+ {
+ initialize(N);
+ set(value);
+ }
+
+ Vector(Subscript N, const T* v) : v_(0), vm1_(0), n_(0)
+ {
+ initialize(N);
+ copy(v);
+ }
+
+ Vector(Subscript N, const char *s) : v_(0), vm1_(0), n_(0)
+ {
+ initialize(N);
+ std::istringstream ins(s);
+
+ Subscript i;
+
+ for (i=0; i<N; i++)
+ ins >> v_[i];
+ }
+
+
+ // methods
+ //
+ Vector<T>& newsize(Subscript N)
+ {
+ if (n_ == N) return *this;
+
+ destroy();
+ initialize(N);
+
+ return *this;
+ }
+
+
+ // assignments
+ //
+ Vector<T>& operator=(const Vector<T> &A)
+ {
+ if (v_ == A.v_)
+ return *this;
+
+ if (n_ == A.n_) // no need to re-alloc
+ copy(A.v_);
+
+ else
+ {
+ destroy();
+ initialize(A.n_);
+ copy(A.v_);
+ }
+
+ return *this;
+ }
+
+ Vector<T>& operator=(const T& scalar)
+ {
+ set(scalar);
+ return *this;
+ }
+
+ inline Subscript dim() const
+ {
+ return n_;
+ }
+
+ inline Subscript size() const
+ {
+ return n_;
+ }
+
+
+ inline reference operator()(Subscript i)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i <= n_) ;
+ #endif
+ return vm1_[i];
+ }
+
+ inline const_reference operator() (Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i <= n_) ;
+ #endif
+ return vm1_[i];
+ }
+
+ inline reference operator[](Subscript i)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(0<=i);
+ assert(i < n_) ;
+ #endif
+ return v_[i];
+ }
+
+ inline const_reference operator[](Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(0<=i);
+
+
+
+
+
+
+ assert(i < n_) ;
+ #endif
+ return v_[i];
+ }
+
+
+
+ };
+
+
+ /* *************************** I/O ********************************/
+
+ template <class T>
+ std::ostream& operator<<(std::ostream &s, const Vector<T> &A)
+ {
+ Subscript N=A.dim();
+
+ s << N << endl;
+
+ for (Subscript i=0; i<N; i++)
+ s << A[i] << " " << endl;
+ s << endl;
+
+ return s;
+ }
+
+ template <class T>
+ std::istream & operator>>(std::istream &s, Vector<T> &A)
+ {
+
+ Subscript N;
+
+ s >> N;
+
+ if ( !(N == A.size() ))
+ {
+ A.newsize(N);
+ }
+
+
+ for (Subscript i=0; i<N; i++)
+ s >> A[i];
+
+
+ return s;
+ }
+
+ // *******************[ basic matrix algorithms ]***************************
+
+
+ template <class T>
+ Vector<T> operator+(const Vector<T> &A,
+ const Vector<T> &B)
+ {
+ Subscript N = A.dim();
+
+ assert(N==B.dim());
+
+ Vector<T> tmp(N);
+ Subscript i;
+
+ for (i=0; i<N; i++)
+ tmp[i] = A[i] + B[i];
+
+ return tmp;
+ }
+
+ template <class T>
+ Vector<T> operator-(const Vector<T> &A,
+ const Vector<T> &B)
+ {
+ Subscript N = A.dim();
+
+ assert(N==B.dim());
+
+ Vector<T> tmp(N);
+ Subscript i;
+
+ for (i=0; i<N; i++)
+ tmp[i] = A[i] - B[i];
+
+ return tmp;
+ }
+
+ template <class T>
+ Vector<T> operator*(const Vector<T> &A,
+ const Vector<T> &B)
+ {
+ Subscript N = A.dim();
+
+ assert(N==B.dim());
+
+ Vector<T> tmp(N);
+ Subscript i;
+
+ for (i=0; i<N; i++)
+ tmp[i] = A[i] * B[i];
+
+ return tmp;
+ }
+
+
+ template <class T>
+ T dot_prod(const Vector<T> &A, const Vector<T> &B)
+ {
+ Subscript N = A.dim();
+ assert(N == B.dim());
+
+ Subscript i;
+ T sum = 0;
+
+ for (i=0; i<N; i++)
+ sum += A[i] * B[i];
+
+ return sum;
+ }
+
+ } /* namespace TNT */
+
+ #endif
+ // VEC_H
Index: tnt/vec.h.patch
===================================================================
RCS file: vec.h.patch
diff -N vec.h.patch
*** /dev/null Tue May 5 14:32:27 1998
--- vec.h.patch Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,15 ----
+ *** vec.h.~1~ Tue Aug 7 15:36:22 2001
+ --- vec.h Tue Aug 7 19:42:18 2001
+ ***************
+ *** 58,62 ****
+ protected:
+ T* v_;
+ ! T* vm1_; // pointer adjustment for optimzied 1-offset indexing
+ Subscript n_;
+
+ --- 58,62 ----
+ protected:
+ T* v_;
+ ! T* vm1_; // pointer adjustment for optimized 1-offset indexing
+ Subscript n_;
+
Index: tnt/vecadaptor.h
===================================================================
RCS file: vecadaptor.h
diff -N vecadaptor.h
*** /dev/null Tue May 5 14:32:27 1998
--- vecadaptor.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,289 ----
+ /*
+ *
+ * Template Numerical Toolkit (TNT): Linear Algebra Module
+ *
+ * Mathematical and Computational Sciences Division
+ * National Institute of Technology,
+ * Gaithersburg, MD USA
+ *
+ *
+ * This software was developed at the National Institute of Standards and
+ * Technology (NIST) by employees of the Federal Government in the course
+ * of their official duties. Pursuant to title 17 Section 105 of the
+ * United States Code, this software is not subject to copyright protection
+ * and is in the public domain. The Template Numerical Toolkit (TNT) is
+ * an experimental system. NIST assumes no responsibility whatsoever for
+ * its use by other parties, and makes no guarantees, expressed or implied,
+ * about its quality, reliability, or any other characteristic.
+ *
+ * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ * see http://math.nist.gov/tnt for latest updates.
+ *
+ */
+
+
+
+ #ifndef VECADAPTOR_H
+ #define VECADAPTOR_H
+
+ #include <cstdlib>
+ #include <iostream>
+ #include <sstream>
+ #include <cassert>
+
+ #include "tnt/subscript.h"
+
+ #ifdef TNT_USE_REGIONS
+ #include "tnt/region1d.h"
+ #endif
+
+ namespace TNT
+ {
+
+ // see "tntreq.h" for TNT requirements for underlying vector
+ // class. This need NOT be the STL vector<> class, but a subset
+ // that provides minimal services.
+ //
+ // This is a container adaptor that provides the following services.
+ //
+ // o) adds 1-offset operator() access ([] is always 0 offset)
+ // o) adds TNT_BOUNDS_CHECK to () and []
+ // o) adds initialization from strings, e.g. "1.0 2.0 3.0";
+ // o) adds newsize(N) function (does not preserve previous values)
+ // o) adds dim() and dim(1)
+ // o) adds free() function to release memory used by vector
+ // o) adds regions, e.g. A(Index(1,10)) = ...
+ // o) add getVector() method to return adapted container
+ // o) adds simple I/O for ostreams
+
+ template <class BBVec>
+ class Vector_Adaptor
+ {
+
+ public:
+ typedef typename BBVec::value_type T;
+ typedef T value_type;
+ typedef T element_type;
+ typedef T* pointer;
+ typedef T* iterator;
+ typedef T& reference;
+ typedef const T* const_iterator;
+ typedef const T& const_reference;
+
+ Subscript lbound() const { return 1; }
+
+ protected:
+ BBVec v_;
+ T* vm1_;
+
+ public:
+
+ Subscript size() const { return v_.size(); }
+
+ // These were removed so that the ANSI C++ valarray class
+ // would work as a possible storage container.
+ //
+ //
+ //iterator begin() { return v_.begin();}
+ //iterator begin() { return &v_[0];}
+ //
+ //iterator end() { return v_.end(); }
+ //iterator end() { return &v_[0] + v_.size(); }
+ //
+ //const_iterator begin() const { return v_.begin();}
+ //const_iterator begin() const { return &v_[0];}
+ //
+ //const_iterator end() const { return v_.end(); }
+ //const_iterator end() const { return &v_[0] + v_.size(); }
+
+ BBVec& getVector() { return v_; }
+ Subscript dim() const { return v_.size(); }
+ Subscript dim(Subscript i)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(i==TNT_BASE_OFFSET);
+ #endif
+ return (i==TNT_BASE_OFFSET ? v_.size() : 0 );
+ }
+ Vector_Adaptor() : v_() {};
+ Vector_Adaptor(const Vector_Adaptor<BBVec> &A) : v_(A.v_)
+ {
+ vm1_ = ( v_.size() > 0 ? &(v_[0]) -1 : NULL);
+
+ }
+
+ Vector_Adaptor(Subscript N, const char *s) : v_(N)
+ {
+ istringstream ins(s);
+ for (Subscript i=0; i<N; i++)
+ ins >> v_[i] ;
+
+ vm1_ = ( v_.size() > 0 ? &(v_[0]) -1 : NULL);
+ };
+
+ Vector_Adaptor(Subscript N, const T& value = T()) : v_(N)
+ {
+ for (Subscript i=0; i<N; i++)
+ v_[i] = value;
+
+ vm1_ = ( v_.size() > 0 ? &(v_[0]) -1 : NULL);
+ }
+
+ Vector_Adaptor(Subscript N, const T* values) : v_(N)
+ {
+ for (Subscript i=0; i<N; i++)
+ v_[i] = values[i];
+ vm1_ = ( v_.size() > 0 ? &(v_[0]) -1 : NULL);
+ }
+ Vector_Adaptor(const BBVec & A) : v_(A)
+ {
+ vm1_ = ( v_.size() > 0 ? &(v_[0]) -1 : NULL);
+ }
+
+ // NOTE: this assumes that BBVec(0) constructor creates an
+ // null vector that does not take up space... It would be
+ // great to require that BBVec have a corresponding free()
+ // function, but in particular STL vectors do not.
+ //
+ Vector_Adaptor<BBVec>& free()
+ {
+ return *this = Vector_Adaptor<BBVec>(0);
+ }
+
+ Vector_Adaptor<BBVec>& operator=(const Vector_Adaptor<BBVec> &A)
+ {
+ v_ = A.v_ ;
+ vm1_ = ( v_.size() > 0 ? &(v_[0]) -1 : NULL);
+ return *this;
+ }
+
+ Vector_Adaptor<BBVec>& newsize(Subscript N)
+ {
+ // NOTE: this is not as efficient as it could be
+ // but to retain compatiblity with STL interface
+ // we cannot assume underlying implementation
+ // has a newsize() function.
+
+ return *this = Vector_Adaptor<BBVec>(N);
+
+ }
+
+ Vector_Adaptor<BBVec>& operator=(const T &a)
+ {
+ Subscript i;
+ Subscript N = v_.size();
+ for (i=0; i<N; i++)
+ v_[i] = a;
+
+ return *this;
+ }
+
+ Vector_Adaptor<BBVec>& resize(Subscript N)
+ {
+ if (N == size()) return *this;
+
+ Vector_Adaptor<BBVec> tmp(N);
+ Subscript n = (N < size() ? N : size()); // min(N, size());
+ Subscript i;
+
+ for (i=0; i<n; i++)
+ tmp[i] = v_[i];
+
+
+ return (*this = tmp);
+
+ }
+
+
+ reference operator()(Subscript i)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i<=dim());
+ #endif
+ return vm1_[i];
+ }
+
+ const_reference operator()(Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(1<=i);
+ assert(i<=dim());
+ #endif
+ return vm1_[i];
+ }
+
+ reference operator[](Subscript i)
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(0<=i);
+ assert(i<dim());
+ #endif
+ return v_[i];
+ }
+
+ const_reference operator[](Subscript i) const
+ {
+ #ifdef TNT_BOUNDS_CHECK
+ assert(0<=i);
+ assert(i<dim());
+ #endif
+ return v_[i];
+ }
+
+
+ #ifdef TNT_USE_REGIONS
+ // "index-aware" features, all of these are 1-based offsets
+
+ typedef Region1D<Vector_Adaptor<BBVec> > Region;
+
+ typedef const_Region1D< Vector_Adaptor<BBVec> > const_Region;
+
+ Region operator()(const Index1D &I)
+ { return Region(*this, I); }
+
+ Region operator()(const Subscript i1, Subscript i2)
+ { return Region(*this, i1, i2); }
+
+ const_Region operator()(const Index1D &I) const
+ { return const_Region(*this, I); }
+
+ const_Region operator()(const Subscript i1, Subscript i2) const
+ { return const_Region(*this, i1, i2); }
+ #endif
+ // TNT_USE_REGIONS
+
+
+ };
+
+ #include <iostream>
+
+ template <class BBVec>
+ std::ostream& operator<<(std::ostream &s, const Vector_Adaptor<BBVec> &A)
+ {
+ Subscript M=A.size();
+
+ s << M << endl;
+ for (Subscript i=1; i<=M; i++)
+ s << A(i) << endl;
+ return s;
+ }
+
+ template <class BBVec>
+ std::istream& operator>>(std::istream &s, Vector_Adaptor<BBVec> &A)
+ {
+ Subscript N;
+
+ s >> N;
+
+ A.resize(N);
+
+ for (Subscript i=1; i<=N; i++)
+ s >> A(i);
+
+ return s;
+ }
+
+ } // namespace TNT
+
+ #endif
Index: tnt/version.h
===================================================================
RCS file: version.h
diff -N version.h
*** /dev/null Tue May 5 14:32:27 1998
--- version.h Thu Aug 16 13:02:56 2001
***************
*** 0 ****
--- 1,20 ----
+ // Template Numerical Toolkit (TNT) for Linear Algebra
+ //
+ // BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
+ // Please see http://math.nist.gov/tnt for updates
+ //
+ // R. Pozo
+ // Mathematical and Computational Sciences Division
+ // National Institute of Standards and Technology
+
+
+ #ifndef TNT_VERSION_H
+ #define TNT_VERSION_H
+
+
+ #define TNT_MAJOR_VERSION '0'
+ #define TNT_MINOR_VERSION '9'
+ #define TNT_SUBMINOR_VERSION '4'
+ #define TNT_VERSION_STRING "0.9.4"
+
+ #endif
More information about the pooma-dev
mailing list