ajout de copies des includes et des sources externes utilisées par le projet Herezh

This commit is contained in:
Gérard Rio 2023-06-15 17:43:43 +02:00
parent 25ae9a89ae
commit 0107e6a9ac
211 changed files with 46684 additions and 0 deletions

Binary file not shown.

File diff suppressed because one or more lines are too long

View file

@ -0,0 +1,876 @@
#!/bin/sh
# This is a shell archive (produced by shar 3.49)
# To extract the files from this archive, save it to a file, remove
# everything above the "!/bin/sh" line above, and type "sh file_name".
#
# made 01/05/1995 18:23 UTC by pozo@salsa
# Source directory /tmp_mnt/home/fs3b/pozo/projects/SparseLib++/1.3/iml/include
#
# existing files will NOT be overwritten unless -c is specified
#
# This shar contains:
# length mode name
# ------ ---------- ------------------------------------------
# 2037 -r--r--r-- bicg.h
# 2227 -r--r--r-- bicgstab.h
# 1711 -r--r--r-- cg.h
# 1978 -r--r--r-- cgs.h
# 2092 -r--r--r-- cheby.h
# 3400 -r--r--r-- gmres.h
# 1379 -r--r--r-- ir.h
# 4089 -r--r--r-- qmr.h
#
# ============= bicg.h ==============
if test -f 'bicg.h' -a X"$1" != X"-c"; then
echo 'x - skipping bicg.h (File already exists)'
else
echo 'x - extracting bicg.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'bicg.h' &&
//*****************************************************************
// Iterative template routine -- BiCG
//
// BiCG solves the unsymmetric linear system Ax = b
// using the Preconditioned BiConjugate Gradient method
//
// BiCG follows the algorithm described on p. 22 of the
// SIAM Templates book.
//
// The return value indicates convergence within max_iter (input)
// iterations (0), or no convergence within max_iter iterations (1).
//
// Upon successful return, output arguments have the following values:
//
// x -- approximate solution to Ax = b
// max_iter -- the number of iterations performed before the
// tolerance was reached
// tol -- the residual after the final iteration
//
//*****************************************************************
X
X
template < class Matrix, class Vector, class Preconditioner, class Real >
int
BiCG(const Matrix &A, Vector &x, const Vector &b,
X const Preconditioner &M, int &max_iter, Real &tol)
{
X Real resid;
X Vector rho_1(1), rho_2(1), alpha(1), beta(1);
X Vector z, ztilde, p, ptilde, q, qtilde;
X
X Real normb = norm(b);
X Vector r = b - A * x;
X Vector rtilde = r;
X
X if (normb == 0.0)
X normb = 1;
X
X if ((resid = norm(r) / normb) <= tol) {
X tol = resid;
X max_iter = 0;
X return 0;
X }
X
X for (int i = 1; i <= max_iter; i++) {
X z = M.solve(r);
X ztilde = M.trans_solve(rtilde);
X rho_1(0) = dot(z, rtilde);
X if (rho_1(0) == 0) {
X tol = norm(r) / normb;
X max_iter = i;
X return 2;
X }
X if (i == 1) {
X p = z;
X ptilde = ztilde;
X } else {
X beta(0) = rho_1(0) / rho_2(0);
X p = z + beta(0) * p;
X ptilde = ztilde + beta(0) * ptilde;
X }
X q = A * p;
X qtilde = A.trans_mult(ptilde);
X alpha(0) = rho_1(0) / dot(ptilde, q);
X x += alpha(0) * p;
X r -= alpha(0) * q;
X rtilde -= alpha(0) * qtilde;
X
X rho_2(0) = rho_1(0);
X if ((resid = norm(r) / normb) < tol) {
X tol = resid;
X max_iter = i;
X return 0;
X }
X }
X
X tol = resid;
X return 1;
}
X
SHAR_EOF
chmod 0444 bicg.h ||
echo 'restore of bicg.h failed'
Wc_c="`wc -c < 'bicg.h'`"
test 2037 -eq "$Wc_c" ||
echo 'bicg.h: original size 2037, current size' "$Wc_c"
fi
# ============= bicgstab.h ==============
if test -f 'bicgstab.h' -a X"$1" != X"-c"; then
echo 'x - skipping bicgstab.h (File already exists)'
else
echo 'x - extracting bicgstab.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'bicgstab.h' &&
//*****************************************************************
// Iterative template routine -- BiCGSTAB
//
// BiCGSTAB solves the unsymmetric linear system Ax = b
// using the Preconditioned BiConjugate Gradient Stabilized method
//
// BiCGSTAB follows the algorithm described on p. 27 of the
// SIAM Templates book.
//
// The return value indicates convergence within max_iter (input)
// iterations (0), or no convergence within max_iter iterations (1).
//
// Upon successful return, output arguments have the following values:
//
// x -- approximate solution to Ax = b
// max_iter -- the number of iterations performed before the
// tolerance was reached
// tol -- the residual after the final iteration
//
//*****************************************************************
X
X
template < class Matrix, class Vector, class Preconditioner, class Real >
int
BiCGSTAB(const Matrix &A, Vector &x, const Vector &b,
X const Preconditioner &M, int &max_iter, Real &tol)
{
X Real resid;
X Vector rho_1(1), rho_2(1), alpha(1), beta(1), omega(1);
X Vector p, phat, s, shat, t, v;
X
X Real normb = norm(b);
X Vector r = b - A * x;
X Vector rtilde = r;
X
X if (normb == 0.0)
X normb = 1;
X
X if ((resid = norm(r) / normb) <= tol) {
X tol = resid;
X max_iter = 0;
X return 0;
X }
X
X for (int i = 1; i <= max_iter; i++) {
X rho_1(0) = dot(rtilde, r);
X if (rho_1(0) == 0) {
X tol = norm(r) / normb;
X return 2;
X }
X if (i == 1)
X p = r;
X else {
X beta(0) = (rho_1(0)/rho_2(0)) * (alpha(0)/omega(0));
X p = r + beta(0) * (p - omega(0) * v);
X }
X phat = M.solve(p);
X v = A * phat;
X alpha(0) = rho_1(0) / dot(rtilde, v);
X s = r - alpha(0) * v;
X if ((resid = norm(s)/normb) < tol) {
X x += alpha(0) * phat;
X tol = resid;
X return 0;
X }
X shat = M.solve(s);
X t = A * shat;
X omega = dot(t,s) / dot(t,t);
X x += alpha(0) * phat + omega(0) * shat;
X r = s - omega(0) * t;
X
X rho_2(0) = rho_1(0);
X if ((resid = norm(r) / normb) < tol) {
X tol = resid;
X max_iter = i;
X return 0;
X }
X if (omega(0) == 0) {
X tol = norm(r) / normb;
X return 3;
X }
X }
X
X tol = resid;
X return 1;
}
SHAR_EOF
chmod 0444 bicgstab.h ||
echo 'restore of bicgstab.h failed'
Wc_c="`wc -c < 'bicgstab.h'`"
test 2227 -eq "$Wc_c" ||
echo 'bicgstab.h: original size 2227, current size' "$Wc_c"
fi
# ============= cg.h ==============
if test -f 'cg.h' -a X"$1" != X"-c"; then
echo 'x - skipping cg.h (File already exists)'
else
echo 'x - extracting cg.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'cg.h' &&
//*****************************************************************
// Iterative template routine -- CG
//
// CG solves the symmetric positive definite linear
// system Ax=b using the Conjugate Gradient method.
//
// CG follows the algorithm described on p. 15 in the
// SIAM Templates book.
//
// The return value indicates convergence within max_iter (input)
// iterations (0), or no convergence within max_iter iterations (1).
//
// Upon successful return, output arguments have the following values:
//
// x -- approximate solution to Ax = b
// max_iter -- the number of iterations performed before the
// tolerance was reached
// tol -- the residual after the final iteration
//
//*****************************************************************
X
template < class Matrix, class Vector, class Preconditioner, class Real >
int
CG(const Matrix &A, Vector &x, const Vector &b,
X const Preconditioner &M, int &max_iter, Real &tol)
{
X Real resid;
X Vector p, z, q;
X Vector alpha(1), beta(1), rho(1), rho_1(1);
X
X Real normb = norm(b);
X Vector r = b - A*x;
X
X if (normb == 0.0)
X normb = 1;
X
X if ((resid = norm(r) / normb) <= tol) {
X tol = resid;
X max_iter = 0;
X return 0;
X }
X
X for (int i = 1; i <= max_iter; i++) {
X z = M.solve(r);
X rho(0) = dot(r, z);
X
X if (i == 1)
X p = z;
X else {
X beta(0) = rho(0) / rho_1(0);
X p = z + beta(0) * p;
X }
X
X q = A*p;
X alpha(0) = rho(0) / dot(p, q);
X
X x += alpha(0) * p;
X r -= alpha(0) * q;
X
X if ((resid = norm(r) / normb) <= tol) {
X tol = resid;
X max_iter = i;
X return 0;
X }
X
X rho_1(0) = rho(0);
X }
X
X tol = resid;
X return 1;
}
X
SHAR_EOF
chmod 0444 cg.h ||
echo 'restore of cg.h failed'
Wc_c="`wc -c < 'cg.h'`"
test 1711 -eq "$Wc_c" ||
echo 'cg.h: original size 1711, current size' "$Wc_c"
fi
# ============= cgs.h ==============
if test -f 'cgs.h' -a X"$1" != X"-c"; then
echo 'x - skipping cgs.h (File already exists)'
else
echo 'x - extracting cgs.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'cgs.h' &&
//*****************************************************************
// Iterative template routine -- CGS
//
// CGS solves the unsymmetric linear system Ax = b
// using the Conjugate Gradient Squared method
//
// CGS follows the algorithm described on p. 26 of the
// SIAM Templates book.
//
// The return value indicates convergence within max_iter (input)
// iterations (0), or no convergence within max_iter iterations (1).
//
// Upon successful return, output arguments have the following values:
//
// x -- approximate solution to Ax = b
// max_iter -- the number of iterations performed before the
// tolerance was reached
// tol -- the residual after the final iteration
//
//*****************************************************************
X
template < class Matrix, class Vector, class Preconditioner, class Real >
int
CGS(const Matrix &A, Vector &x, const Vector &b,
X const Preconditioner &M, int &max_iter, Real &tol)
{
X Real resid;
X Vector rho_1(1), rho_2(1), alpha(1), beta(1);
X Vector p, phat, q, qhat, vhat, u, uhat;
X
X Real normb = norm(b);
X Vector r = b - A*x;
X Vector rtilde = r;
X
X if (normb == 0.0)
X normb = 1;
X
X if ((resid = norm(r) / normb) <= tol) {
X tol = resid;
X max_iter = 0;
X return 0;
X }
X
X for (int i = 1; i <= max_iter; i++) {
X rho_1(0) = dot(rtilde, r);
X if (rho_1(0) == 0) {
X tol = norm(r) / normb;
X return 2;
X }
X if (i == 1) {
X u = r;
X p = u;
X } else {
X beta(0) = rho_1(0) / rho_2(0);
X u = r + beta(0) * q;
X p = u + beta(0) * (q + beta(0) * p);
X }
X phat = M.solve(p);
X vhat = A*phat;
X alpha(0) = rho_1(0) / dot(rtilde, vhat);
X q = u - alpha(0) * vhat;
X uhat = M.solve(u + q);
X x += alpha(0) * uhat;
X qhat = A * uhat;
X r -= alpha(0) * qhat;
X rho_2(0) = rho_1(0);
X if ((resid = norm(r) / normb) < tol) {
X tol = resid;
X max_iter = i;
X return 0;
X }
X }
X
X tol = resid;
X return 1;
}
X
SHAR_EOF
chmod 0444 cgs.h ||
echo 'restore of cgs.h failed'
Wc_c="`wc -c < 'cgs.h'`"
test 1978 -eq "$Wc_c" ||
echo 'cgs.h: original size 1978, current size' "$Wc_c"
fi
# ============= cheby.h ==============
if test -f 'cheby.h' -a X"$1" != X"-c"; then
echo 'x - skipping cheby.h (File already exists)'
else
echo 'x - extracting cheby.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'cheby.h' &&
//*****************************************************************
// Iterative template routine -- CHEBY
//
// CHEBY solves the symmetric positive definite linear
// system Ax = b using the Preconditioned Chebyshev Method
//
// CHEBY follows the algorithm described on p. 30 of the
// SIAM Templates book.
//
// The return value indicates convergence within max_iter (input)
// iterations (0), or no convergence within max_iter iterations (1).
//
// Upon successful return, output arguments have the following values:
//
// x -- approximate solution to Ax = b
// max_iter -- the number of iterations performed before the
// tolerance was reached
// tol -- the residual after the final iteration
//
//*****************************************************************
X
X
template < class Matrix, class Vector, class Preconditioner, class Real,
X class Type >
int
CHEBY(const Matrix &A, Vector &x, const Vector &b,
X const Preconditioner &M, int &max_iter, Real &tol,
X Type eigmin, Type eigmax)
{
X Real resid;
X Type alpha, beta, c, d;
X Vector p, q, z;
X
X Real normb = norm(b);
X Vector r = b - A * x;
X
X if (normb == 0.0)
X normb = 1;
X
X if ((resid = norm(r) / normb) <= tol) {
X tol = resid;
X max_iter = 0;
X return 0;
X }
X
X c = (eigmax - eigmin) / 2.0;
X d = (eigmax + eigmin) / 2.0;
X
X for (int i = 1; i <= max_iter; i++) {
X z = M.solve(r); // apply preconditioner
X
X if (i == 1) {
X p = z;
X alpha = 2.0 / d;
X } else {
X beta = c * alpha / 2.0; // calculate new beta
X beta = beta * beta;
X alpha = 1.0 / (d - beta); // calculate new alpha
X p = z + beta * p; // update search direction
X }
X
X q = A * p;
X x += alpha * p; // update approximation vector
X r -= alpha * q; // compute residual
X
X if ((resid = norm(r) / normb) <= tol) {
X tol = resid;
X max_iter = i;
X return 0; // convergence
X }
X }
X
X tol = resid;
X return 1; // no convergence
}
SHAR_EOF
chmod 0444 cheby.h ||
echo 'restore of cheby.h failed'
Wc_c="`wc -c < 'cheby.h'`"
test 2092 -eq "$Wc_c" ||
echo 'cheby.h: original size 2092, current size' "$Wc_c"
fi
# ============= gmres.h ==============
if test -f 'gmres.h' -a X"$1" != X"-c"; then
echo 'x - skipping gmres.h (File already exists)'
else
echo 'x - extracting gmres.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'gmres.h' &&
//*****************************************************************
// Iterative template routine -- GMRES
//
// GMRES solves the unsymmetric linear system Ax = b using the
// Generalized Minimum Residual method
//
// GMRES follows the algorithm described on p. 20 of the
// SIAM Templates book.
//
// The return value indicates convergence within max_iter (input)
// iterations (0), or no convergence within max_iter iterations (1).
//
// Upon successful return, output arguments have the following values:
//
// x -- approximate solution to Ax = b
// max_iter -- the number of iterations performed before the
// tolerance was reached
// tol -- the residual after the final iteration
//
//*****************************************************************
X
X
template < class Matrix, class Vector >
void
Update(Vector &x, int k, Matrix &h, Vector &s, Vector v[])
{
X Vector y(s);
X
X // Backsolve:
X for (int i = k; i >= 0; i--) {
X y(i) /= h(i,i);
X for (int j = i - 1; j >= 0; j--)
X y(j) -= h(j,i) * y(i);
X }
X
X for (int j = 0; j <= k; j++)
X x += v[j] * y(j);
}
X
X
template < class Real >
Real
abs(Real x)
{
X return (x > 0 ? x : -x);
}
X
X
template < class Operator, class Vector, class Preconditioner,
X class Matrix, class Real >
int
GMRES(const Operator &A, Vector &x, const Vector &b,
X const Preconditioner &M, Matrix &H, int &m, int &max_iter,
X Real &tol)
{
X Real resid;
X int i, j = 1, k;
X Vector s(m+1), cs(m+1), sn(m+1), w;
X
X Real normb = norm(M.solve(b));
X Vector r = M.solve(b - A * x);
X Real beta = norm(r);
X
X if (normb == 0.0)
X normb = 1;
X
X if ((resid = norm(r) / normb) <= tol) {
X tol = resid;
X max_iter = 0;
X return 0;
X }
X
X Vector *v = new Vector[m+1];
X
X while (j <= max_iter) {
X v[0] = r * (1.0 / beta); // ??? r / beta
X s = 0.0;
X s(0) = beta;
X
X for (i = 0; i < m && j <= max_iter; i++, j++) {
X w = M.solve(A * v[i]);
X for (k = 0; k <= i; k++) {
X H(k, i) = dot(w, v[k]);
X w -= H(k, i) * v[k];
X }
X H(i+1, i) = norm(w);
X v[i+1] = w * (1.0 / H(i+1, i)); // ??? w / H(i+1, i)
X
X for (k = 0; k < i; k++)
X ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
X
X GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
X ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
X ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
X
X if ((resid = abs(s(i+1)) / normb) < tol) {
X Update(x, i, H, s, v);
X tol = resid;
X max_iter = j;
X delete [] v;
X return 0;
X }
X }
X Update(x, m - 1, H, s, v);
X r = M.solve(b - A * x);
X beta = norm(r);
X if ((resid = beta / normb) < tol) {
X tol = resid;
X max_iter = j;
X delete [] v;
X return 0;
X }
X }
X
X tol = resid;
X delete [] v;
X return 1;
}
X
X
#include <math.h>
X
X
template<class Real>
void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
{
X if (dy == 0.0) {
X cs = 1.0;
X sn = 0.0;
X } else if (abs(dy) > abs(dx)) {
X Real temp = dx / dy;
X sn = 1.0 / sqrt( 1.0 + temp*temp );
X cs = temp * sn;
X } else {
X Real temp = dy / dx;
X cs = 1.0 / sqrt( 1.0 + temp*temp );
X sn = temp * cs;
X }
}
X
X
template<class Real>
void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
{
X Real temp = cs * dx + sn * dy;
X dy = -sn * dx + cs * dy;
X dx = temp;
}
X
SHAR_EOF
chmod 0444 gmres.h ||
echo 'restore of gmres.h failed'
Wc_c="`wc -c < 'gmres.h'`"
test 3400 -eq "$Wc_c" ||
echo 'gmres.h: original size 3400, current size' "$Wc_c"
fi
# ============= ir.h ==============
if test -f 'ir.h' -a X"$1" != X"-c"; then
echo 'x - skipping ir.h (File already exists)'
else
echo 'x - extracting ir.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'ir.h' &&
//*****************************************************************
// Iterative template routine -- Preconditioned Richardson
//
// IR solves the unsymmetric linear system Ax = b using
// Iterative Refinement (preconditioned Richardson iteration).
//
// The return value indicates convergence within max_iter (input)
// iterations (0), or no convergence within max_iter iterations (1).
//
// Upon successful return, output arguments have the following values:
//
// x -- approximate solution to Ax = b
// max_iter -- the number of iterations performed before the
// tolerance was reached
// tol -- the residual after the final iteration
//
//*****************************************************************
X
template < class Matrix, class Vector, class Preconditioner, class Real >
int
IR(const Matrix &A, Vector &x, const Vector &b,
X const Preconditioner &M, int &max_iter, Real &tol)
{
X Real resid;
X Vector z;
X
X Real normb = norm(b);
X Vector r = b - A*x;
X
X if (normb == 0.0)
X normb = 1;
X
X if ((resid = norm(r) / normb) <= tol) {
X tol = resid;
X max_iter = 0;
X return 0;
X }
X
X for (int i = 1; i <= max_iter; i++) {
X z = M.solve(r);
X x += z;
X r = b - A * x;
X
X if ((resid = norm(r) / normb) <= tol) {
X tol = resid;
X max_iter = i;
X return 0;
X }
X }
X
X tol = resid;
X return 1;
}
X
X
X
SHAR_EOF
chmod 0444 ir.h ||
echo 'restore of ir.h failed'
Wc_c="`wc -c < 'ir.h'`"
test 1379 -eq "$Wc_c" ||
echo 'ir.h: original size 1379, current size' "$Wc_c"
fi
# ============= qmr.h ==============
if test -f 'qmr.h' -a X"$1" != X"-c"; then
echo 'x - skipping qmr.h (File already exists)'
else
echo 'x - extracting qmr.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'qmr.h' &&
//*****************************************************************
// Iterative template routine -- QMR
//
// QMR.h solves the unsymmetric linear system Ax = b using the
// Quasi-Minimal Residual method following the algorithm as described
// on p. 24 in the SIAM Templates book.
//
// -------------------------------------------------------------
// return value indicates
// ------------ ---------------------
// 0 convergence within max_iter iterations
// 1 no convergence after max_iter iterations
// breakdown in:
// 2 rho
// 3 beta
// 4 gamma
// 5 delta
// 6 ep
// 7 xi
// -------------------------------------------------------------
//
// Upon successful return, output arguments have the following values:
//
// x -- approximate solution to Ax=b
// max_iter -- the number of iterations performed before the
// tolerance was reached
// tol -- the residual after the final iteration
//
//*****************************************************************
X
X
#include <math.h>
X
template < class Matrix, class Vector, class Preconditioner1,
X class Preconditioner2, class Real >
int
QMR(const Matrix &A, Vector &x, const Vector &b, const Preconditioner1 &M1,
X const Preconditioner2 &M2, int &max_iter, Real &tol)
{
X Real resid;
X
X Vector rho(1), rho_1(1), xi(1), gamma(1), gamma_1(1), theta(1), theta_1(1);
X Vector eta(1), delta(1), ep(1), beta(1);
X
X Vector r, v_tld, y, w_tld, z;
X Vector v, w, y_tld, z_tld;
X Vector p, q, p_tld, d, s;
X
X Real normb = norm(b);
X
X r = b - A * x;
X
X if (normb == 0.0)
X normb = 1;
X
X if ((resid = norm(r) / normb) <= tol) {
X tol = resid;
X max_iter = 0;
X return 0;
X }
X
X v_tld = r;
X y = M1.solve(v_tld);
X rho(0) = norm(y);
X
X w_tld = r;
X z = M2.trans_solve(w_tld);
X xi(0) = norm(z);
X
X gamma(0) = 1.0;
X eta(0) = -1.0;
X theta(0) = 0.0;
X
X for (int i = 1; i <= max_iter; i++) {
X
X if (rho(0) == 0.0)
X return 2; // return on breakdown
X
X if (xi(0) == 0.0)
X return 7; // return on breakdown
X
X v = (1. / rho(0)) * v_tld;
X y = (1. / rho(0)) * y;
X
X w = (1. / xi(0)) * w_tld;
X z = (1. / xi(0)) * z;
X
X delta(0) = dot(z, y);
X if (delta(0) == 0.0)
X return 5; // return on breakdown
X
X y_tld = M2.solve(y); // apply preconditioners
X z_tld = M1.trans_solve(z);
X
X if (i > 1) {
X p = y_tld - (xi(0) * delta(0) / ep(0)) * p;
X q = z_tld - (rho(0) * delta(0) / ep(0)) * q;
X } else {
X p = y_tld;
X q = z_tld;
X }
X
X p_tld = A * p;
X ep(0) = dot(q, p_tld);
X if (ep(0) == 0.0)
X return 6; // return on breakdown
X
X beta(0) = ep(0) / delta(0);
X if (beta(0) == 0.0)
X return 3; // return on breakdown
X
X v_tld = p_tld - beta(0) * v;
X y = M1.solve(v_tld);
X
X rho_1(0) = rho(0);
X rho(0) = norm(y);
X w_tld = A.trans_mult(q) - beta(0) * w;
X z = M2.trans_solve(w_tld);
X
X xi(0) = norm(z);
X
X gamma_1(0) = gamma(0);
X theta_1(0) = theta(0);
X
X theta(0) = rho(0) / (gamma_1(0) * beta(0));
X gamma(0) = 1.0 / sqrt(1.0 + theta(0) * theta(0));
X
X if (gamma(0) == 0.0)
X return 4; // return on breakdown
X
X eta(0) = -eta(0) * rho_1(0) * gamma(0) * gamma(0) /
X (beta(0) * gamma_1(0) * gamma_1(0));
X
X if (i > 1) {
X d = eta(0) * p + (theta_1(0) * theta_1(0) * gamma(0) * gamma(0)) * d;
X s = eta(0) * p_tld + (theta_1(0) * theta_1(0) * gamma(0) * gamma(0)) * s;
X } else {
X d = eta(0) * p;
X s = eta(0) * p_tld;
X }
X
X x += d; // update approximation vector
X r -= s; // compute residual
X
X if ((resid = norm(r) / normb) <= tol) {
X tol = resid;
X max_iter = i;
X return 0;
X }
X }
X
X tol = resid;
X return 1; // no convergence
}
SHAR_EOF
chmod 0444 qmr.h ||
echo 'restore of qmr.h failed'
Wc_c="`wc -c < 'qmr.h'`"
test 4089 -eq "$Wc_c" ||
echo 'qmr.h: original size 4089, current size' "$Wc_c"
fi
exit 0

File diff suppressed because one or more lines are too long

View file

@ -0,0 +1 @@
//***************************************************************** // Iterative template routine -- BiCG // // BiCG solves the unsymmetric linear system Ax = b // using the Preconditioned BiConjugate Gradient method // // BiCG follows the algorithm described on p. 22 of the // SIAM Templates book. // // The return value indicates convergence within max_iter (input) // iterations (0), or no convergence within max_iter iterations (1). // // Upon successful return, output arguments have the following values: // // x -- approximate solution to Ax = b // max_iter -- the number of iterations performed before the // tolerance was reached // tol -- the residual after the final iteration // //***************************************************************** template < class Matrix, class Vector, class Preconditioner, class Real > int BiCG(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol) { Real resid; Vector rho_1(1), rho_2(1), alpha(1), beta(1); Vector z, ztilde, p, ptilde, q, qtilde; Real normb = norm(b); Vector r = b - A * x; Vector rtilde = r; if (normb == 0.0) normb = 1; if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = 0; return 0; } for (int i = 1; i <= max_iter; i++) { z = M.solve(r); ztilde = M.trans_solve(rtilde); rho_1(0) = dot(z, rtilde); if (rho_1(0) == 0) { tol = norm(r) / normb; max_iter = i; return 2; } if (i == 1) { p = z; ptilde = ztilde; } else { beta(0) = rho_1(0) / rho_2(0); p = z + beta(0) * p; ptilde = ztilde + beta(0) * ptilde; } q = A * p; qtilde = A.trans_mult(ptilde); alpha(0) = rho_1(0) / dot(ptilde, q); x += alpha(0) * p; r -= alpha(0) * q; rtilde -= alpha(0) * qtilde; rho_2(0) = rho_1(0); if ((resid = norm(r) / normb) < tol) { tol = resid; max_iter = i; return 0; } } tol = resid; return 1; }

View file

@ -0,0 +1 @@
//***************************************************************** // Iterative template routine -- BiCGSTAB // // BiCGSTAB solves the unsymmetric linear system Ax = b // using the Preconditioned BiConjugate Gradient Stabilized method // // BiCGSTAB follows the algorithm described on p. 27 of the // SIAM Templates book. // // The return value indicates convergence within max_iter (input) // iterations (0), or no convergence within max_iter iterations (1). // // Upon successful return, output arguments have the following values: // // x -- approximate solution to Ax = b // max_iter -- the number of iterations performed before the // tolerance was reached // tol -- the residual after the final iteration // //***************************************************************** template < class Matrix, class Vector, class Preconditioner, class Real > int BiCGSTAB(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol) { Real resid; Vector rho_1(1), rho_2(1), alpha(1), beta(1), omega(1); Vector p, phat, s, shat, t, v; Real normb = norm(b); Vector r = b - A * x; Vector rtilde = r; if (normb == 0.0) normb = 1; if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = 0; return 0; } for (int i = 1; i <= max_iter; i++) { rho_1(0) = dot(rtilde, r); if (rho_1(0) == 0) { tol = norm(r) / normb; return 2; } if (i == 1) p = r; else { beta(0) = (rho_1(0)/rho_2(0)) * (alpha(0)/omega(0)); p = r + beta(0) * (p - omega(0) * v); } phat = M.solve(p); v = A * phat; alpha(0) = rho_1(0) / dot(rtilde, v); s = r - alpha(0) * v; if ((resid = norm(s)/normb) < tol) { x += alpha(0) * phat; tol = resid; max_iter = i; return 0; } shat = M.solve(s); t = A * shat; omega = dot(t,s) / dot(t,t); x += alpha(0) * phat + omega(0) * shat; r = s - omega(0) * t; rho_2(0) = rho_1(0); if ((resid = norm(r) / normb) < tol) { tol = resid; max_iter = i; return 0; } if (omega(0) == 0) { tol = norm(r) / normb; return 3; } } tol = resid; return 1; }

View file

@ -0,0 +1 @@
//***************************************************************** // Iterative template routine -- CG // // CG solves the symmetric positive definite linear // system Ax=b using the Conjugate Gradient method. // // CG follows the algorithm described on p. 15 in the // SIAM Templates book. // // The return value indicates convergence within max_iter (input) // iterations (0), or no convergence within max_iter iterations (1). // // Upon successful return, output arguments have the following values: // // x -- approximate solution to Ax = b // max_iter -- the number of iterations performed before the // tolerance was reached // tol -- the residual after the final iteration // //***************************************************************** template < class Matrix, class Vector, class Preconditioner, class Real > int CG(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol) { Real resid; Vector p, z, q; Vector alpha(1), beta(1), rho(1), rho_1(1); Real normb = norm(b); Vector r = b - A*x; if (normb == 0.0) normb = 1; if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = 0; return 0; } for (int i = 1; i <= max_iter; i++) { z = M.solve(r); rho(0) = dot(r, z); if (i == 1) p = z; else { beta(0) = rho(0) / rho_1(0); p = z + beta(0) * p; } q = A*p; alpha(0) = rho(0) / dot(p, q); x += alpha(0) * p; r -= alpha(0) * q; if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = i; return 0; } rho_1(0) = rho(0); } tol = resid; return 1; }

View file

@ -0,0 +1 @@
//***************************************************************** // Iterative template routine -- CGS // // CGS solves the unsymmetric linear system Ax = b // using the Conjugate Gradient Squared method // // CGS follows the algorithm described on p. 26 of the // SIAM Templates book. // // The return value indicates convergence within max_iter (input) // iterations (0), or no convergence within max_iter iterations (1). // // Upon successful return, output arguments have the following values: // // x -- approximate solution to Ax = b // max_iter -- the number of iterations performed before the // tolerance was reached // tol -- the residual after the final iteration // //***************************************************************** template < class Matrix, class Vector, class Preconditioner, class Real > int CGS(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol) { Real resid; Vector rho_1(1), rho_2(1), alpha(1), beta(1); Vector p, phat, q, qhat, vhat, u, uhat; Real normb = norm(b); Vector r = b - A*x; Vector rtilde = r; if (normb == 0.0) normb = 1; if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = 0; return 0; } for (int i = 1; i <= max_iter; i++) { rho_1(0) = dot(rtilde, r); if (rho_1(0) == 0) { tol = norm(r) / normb; return 2; } if (i == 1) { u = r; p = u; } else { beta(0) = rho_1(0) / rho_2(0); u = r + beta(0) * q; p = u + beta(0) * (q + beta(0) * p); } phat = M.solve(p); vhat = A*phat; alpha(0) = rho_1(0) / dot(rtilde, vhat); q = u - alpha(0) * vhat; uhat = M.solve(u + q); x += alpha(0) * uhat; qhat = A * uhat; r -= alpha(0) * qhat; rho_2(0) = rho_1(0); if ((resid = norm(r) / normb) < tol) { tol = resid; max_iter = i; return 0; } } tol = resid; return 1; }

View file

@ -0,0 +1 @@
//***************************************************************** // Iterative template routine -- CHEBY // // CHEBY solves the symmetric positive definite linear // system Ax = b using the Preconditioned Chebyshev Method // // CHEBY follows the algorithm described on p. 30 of the // SIAM Templates book. // // The return value indicates convergence within max_iter (input) // iterations (0), or no convergence within max_iter iterations (1). // // Upon successful return, output arguments have the following values: // // x -- approximate solution to Ax = b // max_iter -- the number of iterations performed before the // tolerance was reached // tol -- the residual after the final iteration // //***************************************************************** template < class Matrix, class Vector, class Preconditioner, class Real, class Type > int CHEBY(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol, Type eigmin, Type eigmax) { Real resid; Type alpha, beta, c, d; Vector p, q, z; Real normb = norm(b); Vector r = b - A * x; if (normb == 0.0) normb = 1; if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = 0; return 0; } c = (eigmax - eigmin) / 2.0; d = (eigmax + eigmin) / 2.0; for (int i = 1; i <= max_iter; i++) { z = M.solve(r); // apply preconditioner if (i == 1) { p = z; alpha = 2.0 / d; } else { beta = c * alpha / 2.0; // calculate new beta beta = beta * beta; alpha = 1.0 / (d - beta); // calculate new alpha p = z + beta * p; // update search direction } q = A * p; x += alpha * p; // update approximation vector r -= alpha * q; // compute residual if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = i; return 0; // convergence } } tol = resid; return 1; // no convergence }

View file

@ -0,0 +1 @@
//***************************************************************** // Iterative template routine -- GMRES // // GMRES solves the unsymmetric linear system Ax = b using the // Generalized Minimum Residual method // // GMRES follows the algorithm described on p. 20 of the // SIAM Templates book. // // The return value indicates convergence within max_iter (input) // iterations (0), or no convergence within max_iter iterations (1). // // Upon successful return, output arguments have the following values: // // x -- approximate solution to Ax = b // max_iter -- the number of iterations performed before the // tolerance was reached // tol -- the residual after the final iteration // //***************************************************************** template < class Matrix, class Vector > void Update(Vector &x, int k, Matrix &h, Vector &s, Vector v[]) { Vector y(s); // Backsolve: for (int i = k; i >= 0; i--) { y(i) /= h(i,i); for (int j = i - 1; j >= 0; j--) y(j) -= h(j,i) * y(i); } for (int j = 0; j <= k; j++) x += v[j] * y(j); } template < class Real > Real abs(Real x) { return (x > 0 ? x : -x); } #include <math.h> template<class Real> void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn) { if (dy == 0.0) { cs = 1.0; sn = 0.0; } else if (abs(dy) > abs(dx)) { Real temp = dx / dy; sn = 1.0 / sqrt( 1.0 + temp*temp ); cs = temp * sn; } else { Real temp = dy / dx; cs = 1.0 / sqrt( 1.0 + temp*temp ); sn = temp * cs; } } template<class Real> void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn) { Real temp = cs * dx + sn * dy; dy = -sn * dx + cs * dy; dx = temp; } template < class Operator, class Vector, class Preconditioner, class Matrix, class Real > int GMRES(const Operator &A, Vector &x, const Vector &b, const Preconditioner &M, Matrix &H, int &m, int &max_iter, Real &tol) { Real resid; int i, j = 1, k; Vector s(m+1), cs(m+1), sn(m+1), w; Real normb = norm(M.solve(b)); Vector r = M.solve(b - A * x); Real beta = norm(r); if (normb == 0.0) normb = 1; if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = 0; return 0; } Vector *v = new Vector[m+1]; while (j <= max_iter) { v[0] = r * (1.0 / beta); // ??? r / beta s = 0.0; s(0) = beta; for (i = 0; i < m && j <= max_iter; i++, j++) { w = M.solve(A * v[i]); for (k = 0; k <= i; k++) { H(k, i) = dot(w, v[k]); w -= H(k, i) * v[k]; } H(i+1, i) = norm(w); v[i+1] = w * (1.0 / H(i+1, i)); // ??? w / H(i+1, i) for (k = 0; k < i; k++) ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k)); GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i)); ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i)); ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i)); if ((resid = abs(s(i+1)) / normb) < tol) { Update(x, i, H, s, v); tol = resid; max_iter = j; delete [] v; return 0; } } Update(x, i - 1, H, s, v); r = M.solve(b - A * x); beta = norm(r); if ((resid = beta / normb) < tol) { tol = resid; max_iter = j; delete [] v; return 0; } } tol = resid; delete [] v; return 1; }

View file

@ -0,0 +1 @@
//***************************************************************** // Iterative template routine -- Preconditioned Richardson // // IR solves the unsymmetric linear system Ax = b using // Iterative Refinement (preconditioned Richardson iteration). // // The return value indicates convergence within max_iter (input) // iterations (0), or no convergence within max_iter iterations (1). // // Upon successful return, output arguments have the following values: // // x -- approximate solution to Ax = b // max_iter -- the number of iterations performed before the // tolerance was reached // tol -- the residual after the final iteration // //***************************************************************** template < class Matrix, class Vector, class Preconditioner, class Real > int IR(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol) { Real resid; Vector z; Real normb = norm(b); Vector r = b - A*x; if (normb == 0.0) normb = 1; if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = 0; return 0; } for (int i = 1; i <= max_iter; i++) { z = M.solve(r); x += z; r = b - A * x; if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = i; return 0; } } tol = resid; return 1; }

View file

@ -0,0 +1 @@
//***************************************************************** // Iterative template routine -- QMR // // QMR.h solves the unsymmetric linear system Ax = b using the // Quasi-Minimal Residual method following the algorithm as described // on p. 24 in the SIAM Templates book. // // ------------------------------------------------------------- // return value indicates // ------------ --------------------- // 0 convergence within max_iter iterations // 1 no convergence after max_iter iterations // breakdown in: // 2 rho // 3 beta // 4 gamma // 5 delta // 6 ep // 7 xi // ------------------------------------------------------------- // // Upon successful return, output arguments have the following values: // // x -- approximate solution to Ax=b // max_iter -- the number of iterations performed before the // tolerance was reached // tol -- the residual after the final iteration // //***************************************************************** #include <math.h> template < class Matrix, class Vector, class Preconditioner1, class Preconditioner2, class Real > int QMR(const Matrix &A, Vector &x, const Vector &b, const Preconditioner1 &M1, const Preconditioner2 &M2, int &max_iter, Real &tol) { Real resid; Vector rho(1), rho_1(1), xi(1), gamma(1), gamma_1(1), theta(1), theta_1(1); Vector eta(1), delta(1), ep(1), beta(1); Vector r, v_tld, y, w_tld, z; Vector v, w, y_tld, z_tld; Vector p, q, p_tld, d, s; Real normb = norm(b); r = b - A * x; if (normb == 0.0) normb = 1; if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = 0; return 0; } v_tld = r; y = M1.solve(v_tld); rho(0) = norm(y); w_tld = r; z = M2.trans_solve(w_tld); xi(0) = norm(z); gamma(0) = 1.0; eta(0) = -1.0; theta(0) = 0.0; for (int i = 1; i <= max_iter; i++) { if (rho(0) == 0.0) return 2; // return on breakdown if (xi(0) == 0.0) return 7; // return on breakdown v = (1. / rho(0)) * v_tld; y = (1. / rho(0)) * y; w = (1. / xi(0)) * w_tld; z = (1. / xi(0)) * z; delta(0) = dot(z, y); if (delta(0) == 0.0) return 5; // return on breakdown y_tld = M2.solve(y); // apply preconditioners z_tld = M1.trans_solve(z); if (i > 1) { p = y_tld - (xi(0) * delta(0) / ep(0)) * p; q = z_tld - (rho(0) * delta(0) / ep(0)) * q; } else { p = y_tld; q = z_tld; } p_tld = A * p; ep(0) = dot(q, p_tld); if (ep(0) == 0.0) return 6; // return on breakdown beta(0) = ep(0) / delta(0); if (beta(0) == 0.0) return 3; // return on breakdown v_tld = p_tld - beta(0) * v; y = M1.solve(v_tld); rho_1(0) = rho(0); rho(0) = norm(y); w_tld = A.trans_mult(q) - beta(0) * w; z = M2.trans_solve(w_tld); xi(0) = norm(z); gamma_1(0) = gamma(0); theta_1(0) = theta(0); theta(0) = rho(0) / (gamma_1(0) * beta(0)); gamma(0) = 1.0 / sqrt(1.0 + theta(0) * theta(0)); if (gamma(0) == 0.0) return 4; // return on breakdown eta(0) = -eta(0) * rho_1(0) * gamma(0) * gamma(0) / (beta(0) * gamma_1(0) * gamma_1(0)); if (i > 1) { d = eta(0) * p + (theta_1(0) * theta_1(0) * gamma(0) * gamma(0)) * d; s = eta(0) * p_tld + (theta_1(0) * theta_1(0) * gamma(0) * gamma(0)) * s; } else { d = eta(0) * p; s = eta(0) * p_tld; } x += d; // update approximation vector r -= s; // compute residual if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = i; return 0; } } tol = resid; return 1; // no convergence }

Binary file not shown.

Binary file not shown.

File diff suppressed because it is too large Load diff

File diff suppressed because one or more lines are too long

File diff suppressed because it is too large Load diff

File diff suppressed because one or more lines are too long

Binary file not shown.

Binary file not shown.

File diff suppressed because one or more lines are too long

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ MV++ : Numerical Matrix/Vector Classes in C++ (http://math.nist.gov/pozo/mv++.html) MV++ is a small, efficient, set of concrete vector and simple matrix classes for numerical computing. It is not intended as a general vector container class, but rather designed specifically for optimized numerical computations on RISC and pipelined architectures. It is one step above a C/C++ array, altough it supports subvector/matrix expressions. It includes several of the computational kernels found in the Basic Linear Algebra Subprograms (BLAS), such as scalar updates, vector sums, dot products and so on. These classes are the building blocks of larger user-level libraries such as SparseLib++ and Lapack++. 1) More About MV++ ------------------ i) User's Guide and man-style pages are in http://math.nist.gov/pozo/mv++.html or are available via anonymous ftp from math.nist.gov:pub/pozo/docs/mv++.ps.Z. ii) Code examples are in ./testing/*.cc 2a) Installating MV++ library (float, double, and complex) ------------------------------------------ i) cd to root directory where MV++ is installed, <mv++> ii) edit makefile.def to specify your specify your particular C++ compiler iii) type "make mv"; "make" by itself will provide a list of options. 3) Testing MV++ --------------- i) cd to <mv++> ii) "make test" will run a test suites and leave their output in <mv++>/testing/mv_test.out. 4) Using MV++ ---------- i) all MV++ include files are in <mv++>/include. ii) MV++ library is in <mv++>/lib/mvlib.a iii) to extend non-templated MV++ for vector and matrices of user-defined classes, see mvvt.h, mvmt.h in <mv++>/include. These can easily modified with an editor, replacing "$TYPE" with your specific class name. 5) Help! -------- Questions, comments, suggestions, etc. can be sent to pozo@cam.nist.gov. Questions and Answers: --------------------- o) What types does the non-templated MV++ classes support? Float, double, int, and (optionally) complex. o) How do I extend the non-templated MV++ classes to my own data types? Edit mvvt.h in <mv++>/include and change every occurence of "$TYPE" to the name of your numerical class. (Don't use cpp or m4 for this, as you'll need to change "$TYPE" in word segments as well.) One easy way is sed '1,$s/$$TYPE/Large_Real/g' mvvt.h > mvvLR.h Do the same for mvvt.cc in <mv++>/src. For example, an arbitrary-precision numerical class, Large_Real, would have MV++ classes MV_Vector_Large_Real, MV_ColMat_Large_Real, and so on. o) what is the file "empty" in <mv++>/lib? A dummy file; it was put there so our version of shar would include the <mv++>/lib directory. Otherwise, make reports an error when trying to build mvlib.a , because the destination directory does not exist.

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // read and write MV vectors as text files with (at most) one element // per line. #ifndef _IOTEXT_H_ #define _IOTEXT_H_ int readtxtfile_vec(const char *filename, MV_Vector_double *Aptr); int writetxtfile_vec(const char *filename, const MV_Vector_double &A); int readtxtfile_vec(const char *filename, MV_Vector_int *Aptr); int writetxtfile_vec(const char *filename, const MV_Vector_int &A); #endif

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #ifndef _MV_BLAS1_TPL_H_ #define _MV_BLAS1_TPL_H_ template <class TYPE> MV_Vector<TYPE>& operator*=(MV_Vector<TYPE> &x, const TYPE &a) { int N = x.size(); for (int i=0;i<N;i++) x(i) *= a; return x; } template <class TYPE> MV_Vector<TYPE> operator*(const TYPE &a, const MV_Vector<TYPE> &x) { int N = x.size(); MV_Vector<TYPE> result(N); for (int i=0;i<N;i++) result(i) = x(i)*a; return result; } template <class TYPE> MV_Vector<TYPE> operator*(const MV_Vector<TYPE> &x, const TYPE &a) { // This is the other commutative case of vector*scalar. // It should be just defined to be // "return operator*(a,x);" // but some compilers (e.g. Turbo C++ v.3.0) have trouble // determining the proper template match. For the moment, // we'll just duplicate the code in the scalar * MV_Vector // case above. int N = x.size(); MV_Vector<TYPE> result(N); for (int i=0;i<N;i++) result(i) = x(i)*a; return result; } template <class TYPE> MV_Vector<TYPE> operator+(const MV_Vector<TYPE> &x, const MV_Vector<TYPE> &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in +." << endl; exit(1); } MV_Vector<TYPE> result(N); for (int i=0;i<N; i++) result(i) = x(i) + y(i); return result; } template <class TYPE> MV_Vector<TYPE> operator-(const MV_Vector<TYPE> &x, const MV_Vector<TYPE> &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } MV_Vector<TYPE> result(N); for (int i=0;i<N; i++) result(i) = x(i) - y(i); return result; } template <class TYPE> MV_Vector<TYPE>& operator+=(MV_Vector<TYPE> &x, const MV_Vector<TYPE> &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } for (int i=0;i<N; i++) x(i) += y(i); return x; } template <class TYPE> MV_Vector<TYPE>& operator-=(MV_Vector<TYPE> &x, const MV_Vector<TYPE> &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } for (int i=0;i<N; i++) x(i) -= y(i); return x; } // norm and dot product functions for the MV_Vector<> class template <class TYPE> TYPE dot(const MV_Vector<TYPE> &x, const MV_Vector<TYPE> &y) { // Check for compatible dimensions: if (x.size() != y.size()) { cout << "Incompatible dimensions in dot(). " << endl; exit(1); } TYPE temp=0.0; for (int i=0; i<x.size();i++) temp += x(i)*y(i); return temp; } template <class TYPE> TYPE norm(const MV_Vector<TYPE> &x) { TYPE temp = dot(x,x); return sqrt(temp); } #endif // _MV_BLAS1_TPL_H_

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #ifndef _MV_BLAS1_COMPLEX_H_ #define _MV_BLAS1_COMPLEX_H_ MV_Vector_COMPLEX& operator*=(MV_Vector_COMPLEX &x, const COMPLEX &a); MV_Vector_COMPLEX operator*(const COMPLEX &a, const MV_Vector_COMPLEX &x); MV_Vector_COMPLEX operator*(const MV_Vector_COMPLEX &x, const COMPLEX &a); MV_Vector_COMPLEX operator+(const MV_Vector_COMPLEX &x, const MV_Vector_COMPLEX &y); MV_Vector_COMPLEX operator-(const MV_Vector_COMPLEX &x, const MV_Vector_COMPLEX &y); MV_Vector_COMPLEX& operator+=(MV_Vector_COMPLEX &x, const MV_Vector_COMPLEX &y); MV_Vector_COMPLEX& operator-=(MV_Vector_COMPLEX &x, const MV_Vector_COMPLEX &y); COMPLEX dot(const MV_Vector_COMPLEX &x, const MV_Vector_COMPLEX &y); COMPLEX norm(const MV_Vector_COMPLEX &x); #endif // _MV_BLAS1_COMPLEX_H_

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #ifndef _MV_BLAS1_double_H_ #define _MV_BLAS1_double_H_ MV_Vector_double& operator*=(MV_Vector_double &x, const double &a); MV_Vector_double operator*(const double &a, const MV_Vector_double &x); MV_Vector_double operator*(const MV_Vector_double &x, const double &a); MV_Vector_double operator+(const MV_Vector_double &x, const MV_Vector_double &y); MV_Vector_double operator-(const MV_Vector_double &x, const MV_Vector_double &y); MV_Vector_double& operator+=(MV_Vector_double &x, const MV_Vector_double &y); MV_Vector_double& operator-=(MV_Vector_double &x, const MV_Vector_double &y); double dot(const MV_Vector_double &x, const MV_Vector_double &y); double norm(const MV_Vector_double &x); #endif // _MV_BLAS1_double_H_

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #ifndef _MV_BLAS1_float_H_ #define _MV_BLAS1_float_H_ MV_Vector_float& operator*=(MV_Vector_float &x, const float &a); MV_Vector_float operator*(const float &a, const MV_Vector_float &x); MV_Vector_float operator*(const MV_Vector_float &x, const float &a); MV_Vector_float operator+(const MV_Vector_float &x, const MV_Vector_float &y); MV_Vector_float operator-(const MV_Vector_float &x, const MV_Vector_float &y); MV_Vector_float& operator+=(MV_Vector_float &x, const MV_Vector_float &y); MV_Vector_float& operator-=(MV_Vector_float &x, const MV_Vector_float &y); float dot(const MV_Vector_float &x, const MV_Vector_float &y); float norm(const MV_Vector_float &x); #endif // _MV_BLAS1_float_H_

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #ifndef _MV_BLAS1_int_H_ #define _MV_BLAS1_int_H_ MV_Vector_int& operator*=(MV_Vector_int &x, const int &a); MV_Vector_int operator*(const int &a, const MV_Vector_int &x); MV_Vector_int operator*(const MV_Vector_int &x, const int &a); MV_Vector_int operator+(const MV_Vector_int &x, const MV_Vector_int &y); MV_Vector_int operator-(const MV_Vector_int &x, const MV_Vector_int &y); MV_Vector_int& operator+=(MV_Vector_int &x, const MV_Vector_int &y); MV_Vector_int& operator-=(MV_Vector_int &x, const MV_Vector_int &y); int dot(const MV_Vector_int &x, const MV_Vector_int &y); int norm(const MV_Vector_int &x); #endif // _MV_BLAS1_int_H_

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #ifndef _MV_BLAS1_$TYPE_H_ #define _MV_BLAS1_$TYPE_H_ MV_Vector_$TYPE& operator*=(MV_Vector_$TYPE &x, const $TYPE &a); MV_Vector_$TYPE operator*(const $TYPE &a, const MV_Vector_$TYPE &x); MV_Vector_$TYPE operator*(const MV_Vector_$TYPE &x, const $TYPE &a); MV_Vector_$TYPE operator+(const MV_Vector_$TYPE &x, const MV_Vector_$TYPE &y); MV_Vector_$TYPE operator-(const MV_Vector_$TYPE &x, const MV_Vector_$TYPE &y); MV_Vector_$TYPE& operator+=(MV_Vector_$TYPE &x, const MV_Vector_$TYPE &y); MV_Vector_$TYPE& operator-=(MV_Vector_$TYPE &x, const MV_Vector_$TYPE &y); $TYPE dot(const MV_Vector_$TYPE &x, const MV_Vector_$TYPE &y); $TYPE norm(const MV_Vector_$TYPE &x); #endif // _MV_BLAS1_$TYPE_H_

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // mvm.h // // MV Matrix classes for double, float, int, and complex // #ifndef _MV_MATRIX_ALL_H_ #define _MV_MATRIX_ALL_H_ #include "mvmrf.h" #include "mvmd.h" #include "mvmi.h" #include "mvmf.h" #include <complex.h> #include "mvmc.h" #endif // _MV_MATRIX_ALL_H_

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // this is really used as a sort of global constant. The reason // for creating its own type is that so it can be overloaded to perform // a deep or shallow assignement. (Any variable of type MV_Matrix_::ref_type // has only one possible value: one.) #ifndef _MV_MATRIX_REF_ #define _MV_MATRIX_REF_ struct MV_Matrix_ { enum ref_type { ref = 1}; } ; # endif

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // mvv.h // // MV vector classes for double, float, int, and complex // #ifndef _MV_VECTOR_ALL_H_ #define _MV_VECTOR_ALL_H_ #include "mvvrf.h" #include "mvvd.h" #include "mvvi.h" #include "mvvf.h" #include "mvblasf.h" #include "mvblasd.h" #include "mvblasi.h" #include <complex.h> #include "mvvc.h" #include "mvblasc.h" #include "iotext.h" #endif // _MV_VECTOR_ALL_H_

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // mvvind.h MV_Vector Index class #ifndef _MV_VEC_INDEX_H_ #define _MV_VEC_INDEX_H_ // A MV_VecIndex is an ordered pair (start,end) denoting a subvector // region, similar to a Fortran 90 or Matlab colon notation. For example, // // MV_Vector_double A(10), B(20); // MV_VecIndex I(2,4); // // A(I) = B(MV_VecIndex(0,2); // // sets the thrid through fifth elements of A to the first two elements // of B. There is no stride argument, only contiguous regions are allowed. // #include <assert.h> class MV_VecIndex { private: unsigned int start_; unsigned int end_; char all_; // true if this index refers to the complete // vector range. start_ and end_ are ignored. public: MV_VecIndex() : start_(0), end_(0), all_(1) {} MV_VecIndex(unsigned int i1) :start_(i1), end_(i1), all_(0) {} MV_VecIndex(unsigned int i1, unsigned int i2): start_(i1), end_(i2), all_(0) { assert(i1 <= i2); } MV_VecIndex(const MV_VecIndex &s) : start_(s.start_), end_(s.end_), all_(s.all_){} int start() const { return (all_==1) ? 0 : start_;} int end() const { return (all_ ==1) ? 0 : end_;} int length() const { return (all_==1) ? 0 : (end_-start_+1);} int all() const { return all_; } MV_VecIndex& operator=(const MV_VecIndex& I) { start_=I.start_; end_ = I.end_; return *this;} MV_VecIndex operator+(int i) { return MV_VecIndex(start_ +i, end_ +i); } MV_VecIndex& operator+=(int i) { start_ += i; end_ += i; return *this; } MV_VecIndex operator-(int i) { return MV_VecIndex(start_ -i, end_ -i); } MV_VecIndex& operator-=(int i) { start_ -= i; end_ -= i; return *this; } }; #endif // _INDEX_H_

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // this is really used as a sort of global constant. The reason // for creating its own type is that so it can be overloaded to perform // a deep or shallow assignement. (Any variable of type MV_Vector_::ref_type // has only one possible value: one.) #ifndef _MV_VECTOR_REF_ #define _MV_VECTOR_REF_ struct MV_Vector_ { enum ref_type { ref = 1}; } ; # endif

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #ifndef _MV_BLAS1_TYPE_H_ #define _MV_BLAS1_TYPE_H_ #include <math.h> #include <stdlib.h> MV_Vector_TYPE& operator*=(MV_Vector_TYPE &x, const TYPE &a); MV_Vector_TYPE operator*(const TYPE &a, const MV_Vector_TYPE &x); MV_Vector_TYPE operator*(const MV_Vector_TYPE &x, const TYPE &a); MV_Vector_TYPE operator+(const MV_Vector_TYPE &x, const MV_Vector_TYPE &y); MV_Vector_TYPE operator-(const MV_Vector_TYPE &x, const MV_Vector_TYPE &y); MV_Vector_TYPE& operator+=(MV_Vector_TYPE &x, const MV_Vector_TYPE &y); MV_Vector_TYPE& operator-=(MV_Vector_TYPE &x, const MV_Vector_TYPE &y); TYPE dot(const MV_Vector_TYPE &x, const MV_Vector_TYPE &y); TYPE norm(const MV_Vector_TYPE &x); #endif // _MV_BLAS1_TYPE_H_

File diff suppressed because one or more lines are too long

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* Which dense vector/matrix classes to build SparseLib++ from */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #ifndef vector_defs_H #define vector_defs_H #define VECTOR_H "mvv.h" #define VECTOR_double MV_Vector_double #define VECTOR_float MV_Vector_float #define VECTOR_int MV_Vector_int #define VECTOR_ref MV_Vector_::ref #define MATRIX_H "mvm.h" #define MATRIX_double MV_ColMat_double #define MATRIX_float MV_ColMat_float #define MATRIX_int MV_ColMat_int #define MATRIX_ref MV_Matrix_::ref #define VECTOR_COMPLEX MV_Vector_COMPLEX #define MATRIX_COMPLEX MV_ColMat_COMPLEX #endif

View file

@ -0,0 +1 @@
This is a dummy file; it was put there so our version of shar would include the <mv++>/lib directory. Otherwise, make reports an error when trying to build mvlib.a , because the destination directory does not exist.

View file

@ -0,0 +1 @@
################################################################ # # Subdirectories for building # ################################################################ error: @echo "+-----------------------------------------------------------------+" @echo "| |" @echo "| MV++ Matrix/Vector C++ Class Library |" @echo "| |" @echo "| Usage: make all install and test MV++ |" @echo "| |" @echo "| make mvlib install non-templated MV++ library |" @echo "| make mvtest run non-templated test suite |" @echo "| make clean clean *.o and test executables |" @echo "| |" @echo "| Make sure the system-specific makefile.def has been edited |" @echo "| to reflect your system configuration. |" @echo "| |" @echo "+-----------------------------------------------------------------+" all: mvlib mvtest done mvlib: cd ./src; make; mvtest: cd ./testing; make; cd ./testing; mvtest > mvtest.out; @echo " " @echo " +---------------------------------------------------------------+" @echo " | |" @echo " | MV++ test complete. Trace is in <mv++>/testing/mvtest.out |" @echo " | |" @echo " +---------------------------------------------------------------+" @echo " " clean: cd ./src; make clean; cd ./testing; make clean; wipe: cd ./src; make wipe; cd ./testing; make wipe; done: @echo " " @echo " +---------------------------------------------------------------+" @echo " | |" @echo " | MV++ installed. |" @echo " | |" @echo " +---------------------------------------------------------------+" @echo " " @echo " "

View file

@ -0,0 +1 @@
#*************************************************************************** # # M V ++ M A K E F I L E . D E F # # Compiler Specific Section: Edit to reflect your environment #--------------------------------------------------------------------------- # # Macro Significance # # CCC C++ compiler # CCCFLAGS flags to C++ compiler # CC C compiler # CFLAGS flags to C compiler # LDFLAGS libraries # HASRANLIB 't' if your system has ranlib, 'f' otherwise # # # Some defaults are below for Sun C++, and GNU g++. Uncomment the # relevant sections. # # # # Support for optional vector-bounds checking: # # Use -DMV_VECTOR_BOUNDS_CHECK to perform runtime checking of # index bounds. Can be used to create a "debug" and a # a "production" version of MV++ libraries. # # Support for vectors and matrices of type complex: # # Use -DCOMPLEX_OSTREAM if your compiler does not provide a predefine # a way to print complex numbers (e.g. cout << u; ) # # Define COMPLEX to expand to whatever name your compiler uses for # complex numbers (e.g. g++ 2.7.0 uses templated classes # as per ANSI C++) Most other compilers use "complex" or # "Complex". # IBM xlC v. 1.1 # CCC = xlC CC = xlc CFLAGS = -O CCCFLAGS = -+ -O -DCOMPLEX=complex -DMV_VECTOR_BOUNDS_CHECK LDFLAGS = ../lib/libmv.a -lm -lcomplex # g++ v. 2.7.0 # CCC = g++ CC = gcc CFLAGS = -O CCCFLAGS = -DMV_VECTOR_BOUNDS_CHECK -'DCOMPLEX=complex<double>' LDFLAGS = ../lib/libmv.a -lm # g++ v. 2.6.3 # CCC = g++ CC = gcc CFLAGS = -O CCCFLAGS = -Wall -DMV_VECTOR_BOUNDS_CHECK -DCOMPLEX=complex LDFLAGS = ../lib/libmv.a -lm # Sun C++ 4.0.1 # CCC = CC CC = acc CFLAGS = -O CCCFLAGS = +w -DCOMPLEX_OSTREAM -I/usr/local/apps/lang1/SC3.0.1/include/CC_413/ -DCOMPLEX=complex -DMV_VECTOR_BOUNDS_CHECK LDFLAGS = ../lib/libmv.a -lm -lcomplex # ranlib available on this system? 't' or 'f' HASRANLIB = t ################################################################ # # # Implicit rules: do not modfiy after here # # # ################################################################ MV_DIR = .. MV_LIB_DIR = $(MV_DIR)/lib MV_LIB_A = $(MV_LIB_DIR)/libmv.a MV_INCLUDE_DIR = $(MV_DIR)/include EXENAME = -o LIB_EXT = a EXE_EXT = .SUFFIXES: .c .cc .o $(EXE_EXT) .o$(EXE_EXT): $(CCC) $(EXENAME) $@$(EXE_EXT) $@.o \ $(LIBS)$(LDFLAGS) .cc.o: $(CCC) $(CCCFLAGS) $(INCLUDES) -c $< .c.o: $(CC) $(CFLAGS) $(INCLUDES) -c $< .f.o: $(F77) -c $<

View file

@ -0,0 +1 @@
# # Makefile for sparselb src # include ../makefile.def MV_INCLUDE_DIR = ../include INCLUDES = -I$(MV_INCLUDE_DIR) SRCS = \ mvvd.cc \ mvmd.cc \ mvvf.cc \ mvmf.cc \ mvvdio.cc \ mvblasd.cc \ mvblasf.cc \ mvblasi.cc \ mvvi.cc \ mvmi.cc \ mvvc.cc \ mvmc.cc \ mvvcio.cc \ mvblasc.cc OBJS = \ mvvd.o \ mvmd.o \ mvblasd.o \ mvvf.o \ mvmf.o \ mvblasi.o \ mvblasf.o \ mvvi.o \ mvmi.o \ mvvdio.o \ mvvc.o \ mvmc.o \ mvvcio.o \ mvblasc.o all: lib srcs: $(SRCS) lib: $(MV_LIB_A) #$(MV_LIB_A) : $(OBJS) # @ ar rv $(MV_LIB_A) $? # @ case x$(HASRANLIB) in xt ) echo ranlib; ranlib $(MV_LIB_A);; esac $(MV_LIB_A) : $(OBJS) @ ar rv $(MV_LIB_A) $? @ case x$(HASRANLIB) in xt ) echo ranlib; ranlib $(MV_LIB_A);; esac clean: /bin/rm -f $(OBJS) wipe: /bin/rm -f $(OBJS) $(MV_LIB_A)

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #include <math.h> #include <stdlib.h> #include "mvv.h" MV_Vector_COMPLEX& operator*=(MV_Vector_COMPLEX &x, const COMPLEX &a) { int N = x.size(); for (int i=0;i<N;i++) x(i) *= a; return x; } MV_Vector_COMPLEX operator*(const COMPLEX &a, const MV_Vector_COMPLEX &x) { int N = x.size(); MV_Vector_COMPLEX result(N); for (int i=0;i<N;i++) result(i) = x(i)*a; return result; } MV_Vector_COMPLEX operator*(const MV_Vector_COMPLEX &x, const COMPLEX &a) { // This is the other commutative case of vector*scalar. // It should be just defined to be // "return operator*(a,x);" // but some compilers (e.g. Turbo C++ v.3.0) have trouble // determining the proper template match. For the moment, // we'll just duplicate the code in the scalar * vector // case above. int N = x.size(); MV_Vector_COMPLEX result(N); for (int i=0;i<N;i++) result(i) = x(i)*a; return result; } MV_Vector_COMPLEX operator+(const MV_Vector_COMPLEX &x, const MV_Vector_COMPLEX &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in +." << endl; exit(1); } MV_Vector_COMPLEX result(N); for (int i=0;i<N; i++) result(i) = x(i) + y(i); return result; } MV_Vector_COMPLEX operator-(const MV_Vector_COMPLEX &x, const MV_Vector_COMPLEX &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } MV_Vector_COMPLEX result(N); for (int i=0;i<N; i++) result(i) = x(i) - y(i); return result; } MV_Vector_COMPLEX& operator+=(MV_Vector_COMPLEX &x, const MV_Vector_COMPLEX &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } for (int i=0;i<N; i++) x(i) += y(i); return x; } MV_Vector_COMPLEX& operator-=(MV_Vector_COMPLEX &x, const MV_Vector_COMPLEX &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } for (int i=0;i<N; i++) x(i) -= y(i); return x; } // norm and dot product functions for the MV_Vector<> class COMPLEX dot(const MV_Vector_COMPLEX &x, const MV_Vector_COMPLEX &y) { // Check for compatible dimensions: if (x.size() != y.size()) { cout << "Incompatible dimensions in dot(). " << endl; exit(1); } COMPLEX temp = 0; for (int i=0; i<x.size();i++) temp += x(i)*y(i); return temp; } COMPLEX norm(const MV_Vector_COMPLEX &x) { COMPLEX temp = dot(x,x); return (COMPLEX) sqrt(temp); }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #include <math.h> #include <stdlib.h> #include "mvvd.h" MV_Vector_double& operator*=(MV_Vector_double &x, const double &a) { int N = x.size(); for (int i=0;i<N;i++) x(i) *= a; return x; } MV_Vector_double operator*(const double &a, const MV_Vector_double &x) { int N = x.size(); MV_Vector_double result(N); for (int i=0;i<N;i++) result(i) = x(i)*a; return result; } MV_Vector_double operator*(const MV_Vector_double &x, const double &a) { // This is the other commutative case of vector*scalar. // It should be just defined to be // "return operator*(a,x);" // but some compilers (e.g. Turbo C++ v.3.0) have trouble // determining the proper template match. For the moment, // we'll just duplicate the code in the scalar * vector // case above. int N = x.size(); MV_Vector_double result(N); for (int i=0;i<N;i++) result(i) = x(i)*a; return result; } MV_Vector_double operator+(const MV_Vector_double &x, const MV_Vector_double &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in +." << endl; exit(1); } MV_Vector_double result(N); for (int i=0;i<N; i++) result(i) = x(i) + y(i); return result; } MV_Vector_double operator-(const MV_Vector_double &x, const MV_Vector_double &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } MV_Vector_double result(N); for (int i=0;i<N; i++) result(i) = x(i) - y(i); return result; } MV_Vector_double& operator+=(MV_Vector_double &x, const MV_Vector_double &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } for (int i=0;i<N; i++) x(i) += y(i); return x; } MV_Vector_double& operator-=(MV_Vector_double &x, const MV_Vector_double &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } for (int i=0;i<N; i++) x(i) -= y(i); return x; } // norm and dot product functions for the MV_Vector<> class double dot(const MV_Vector_double &x, const MV_Vector_double &y) { // Check for compatible dimensions: if (x.size() != y.size()) { cout << "Incompatible dimensions in dot(). " << endl; exit(1); } double temp = 0; for (int i=0; i<x.size();i++) temp += x(i)*y(i); return temp; } double norm(const MV_Vector_double &x) { double temp = dot(x,x); return (double) sqrt(temp); }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #include <math.h> #include <stdlib.h> #include "mvvf.h" MV_Vector_float& operator*=(MV_Vector_float &x, const float &a) { int N = x.size(); for (int i=0;i<N;i++) x(i) *= a; return x; } MV_Vector_float operator*(const float &a, const MV_Vector_float &x) { int N = x.size(); MV_Vector_float result(N); for (int i=0;i<N;i++) result(i) = x(i)*a; return result; } MV_Vector_float operator*(const MV_Vector_float &x, const float &a) { // This is the other commutative case of vector*scalar. // It should be just defined to be // "return operator*(a,x);" // but some compilers (e.g. Turbo C++ v.3.0) have trouble // determining the proper template match. For the moment, // we'll just duplicate the code in the scalar * vector // case above. int N = x.size(); MV_Vector_float result(N); for (int i=0;i<N;i++) result(i) = x(i)*a; return result; } MV_Vector_float operator+(const MV_Vector_float &x, const MV_Vector_float &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in +." << endl; exit(1); } MV_Vector_float result(N); for (int i=0;i<N; i++) result(i) = x(i) + y(i); return result; } MV_Vector_float operator-(const MV_Vector_float &x, const MV_Vector_float &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } MV_Vector_float result(N); for (int i=0;i<N; i++) result(i) = x(i) - y(i); return result; } MV_Vector_float& operator+=(MV_Vector_float &x, const MV_Vector_float &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } for (int i=0;i<N; i++) x(i) += y(i); return x; } MV_Vector_float& operator-=(MV_Vector_float &x, const MV_Vector_float &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } for (int i=0;i<N; i++) x(i) -= y(i); return x; } // norm and dot product functions for the MV_Vector<> class float dot(const MV_Vector_float &x, const MV_Vector_float &y) { // Check for compatible dimensions: if (x.size() != y.size()) { cout << "Incompatible dimensions in dot(). " << endl; exit(1); } float temp = 0; for (int i=0; i<x.size();i++) temp += x(i)*y(i); return temp; } float norm(const MV_Vector_float &x) { float temp = dot(x,x); return (float) sqrt(temp); }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #include <math.h> #include <stdlib.h> #include "mvvi.h" MV_Vector_int& operator*=(MV_Vector_int &x, const int &a) { int N = x.size(); for (int i=0;i<N;i++) x(i) *= a; return x; } MV_Vector_int operator*(const int &a, const MV_Vector_int &x) { int N = x.size(); MV_Vector_int result(N); for (int i=0;i<N;i++) result(i) = x(i)*a; return result; } MV_Vector_int operator*(const MV_Vector_int &x, const int &a) { // This is the other commutative case of vector*scalar. // It should be just defined to be // "return operator*(a,x);" // but some compilers (e.g. Turbo C++ v.3.0) have trouble // determining the proper template match. For the moment, // we'll just duplicate the code in the scalar * vector // case above. int N = x.size(); MV_Vector_int result(N); for (int i=0;i<N;i++) result(i) = x(i)*a; return result; } MV_Vector_int operator+(const MV_Vector_int &x, const MV_Vector_int &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in +." << endl; exit(1); } MV_Vector_int result(N); for (int i=0;i<N; i++) result(i) = x(i) + y(i); return result; } MV_Vector_int operator-(const MV_Vector_int &x, const MV_Vector_int &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } MV_Vector_int result(N); for (int i=0;i<N; i++) result(i) = x(i) - y(i); return result; } MV_Vector_int& operator+=(MV_Vector_int &x, const MV_Vector_int &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } for (int i=0;i<N; i++) x(i) += y(i); return x; } MV_Vector_int& operator-=(MV_Vector_int &x, const MV_Vector_int &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } for (int i=0;i<N; i++) x(i) -= y(i); return x; } // norm and dot product functions for the MV_Vector<> class int dot(const MV_Vector_int &x, const MV_Vector_int &y) { // Check for compatible dimensions: if (x.size() != y.size()) { cout << "Incompatible dimensions in dot(). " << endl; exit(1); } int temp = 0; for (int i=0; i<x.size();i++) temp += x(i)*y(i); return temp; } int norm(const MV_Vector_int &x) { int temp = dot(x,x); return (int) sqrt(temp); }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #include <math.h> #include <stdlib.h> #include "$INCLUDE" MV_Vector_$TYPE& operator*=(MV_Vector_$TYPE &x, const $TYPE &a) { int N = x.size(); for (int i=0;i<N;i++) x(i) *= a; return x; } MV_Vector_$TYPE operator*(const $TYPE &a, const MV_Vector_$TYPE &x) { int N = x.size(); MV_Vector_$TYPE result(N); for (int i=0;i<N;i++) result(i) = x(i)*a; return result; } MV_Vector_$TYPE operator*(const MV_Vector_$TYPE &x, const $TYPE &a) { // This is the other commutative case of vector*scalar. // It should be just defined to be // "return operator*(a,x);" // but some compilers (e.g. Turbo C++ v.3.0) have trouble // determining the proper template match. For the moment, // we'll just duplicate the code in the scalar * vector // case above. int N = x.size(); MV_Vector_$TYPE result(N); for (int i=0;i<N;i++) result(i) = x(i)*a; return result; } MV_Vector_$TYPE operator+(const MV_Vector_$TYPE &x, const MV_Vector_$TYPE &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in +." << endl; exit(1); } MV_Vector_$TYPE result(N); for (int i=0;i<N; i++) result(i) = x(i) + y(i); return result; } MV_Vector_$TYPE operator-(const MV_Vector_$TYPE &x, const MV_Vector_$TYPE &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } MV_Vector_$TYPE result(N); for (int i=0;i<N; i++) result(i) = x(i) - y(i); return result; } MV_Vector_$TYPE& operator+=(MV_Vector_$TYPE &x, const MV_Vector_$TYPE &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } for (int i=0;i<N; i++) x(i) += y(i); return x; } MV_Vector_$TYPE& operator-=(MV_Vector_$TYPE &x, const MV_Vector_$TYPE &y) { int N = x.size(); if (N != y.size()) { cout << "Incompatible vector lengths in -." << endl; exit(1); } for (int i=0;i<N; i++) x(i) -= y(i); return x; } // norm and dot product functions for the MV_Vector<> class $TYPE dot(const MV_Vector_$TYPE &x, const MV_Vector_$TYPE &y) { // Check for compatible dimensions: if (x.size() != y.size()) { cout << "Incompatible dimensions in dot(). " << endl; exit(1); } $TYPE temp = 0; for (int i=0; i<x.size();i++) temp += x(i)*y(i); return temp; } $TYPE norm(const MV_Vector_$TYPE &x) { $TYPE temp = dot(x,x); return ($TYPE) sqrt(temp); }

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

View file

@ -0,0 +1 @@
// // ostream_complex.cc Add I/O routine for printing out complex // numbers. // #ifdef COMPLEX_OSTREAM #include <complex.h> // AT&T Cfront does not provide for cout << complex ... // ostream& operator<<(ostream &s, COMPLEX z) { s << (double) real(z) << " " << (double) imag(z) ; return s; } #endif

File diff suppressed because one or more lines are too long

View file

@ -0,0 +1 @@
#include <stdio.h> #include "mvvd.h" #include "mvvi.h" #include "iotext.h" int writetxtfile_vec(const char *filename, const MV_Vector_double &A) { FILE *in_file; in_file = fopen(filename, "w"); if (in_file == NULL) { fprintf(stderr,"Cannot open file: %s\n", filename ); exit(1); } int N = A.size(); for (int i=0; i<N; i++) { fprintf(in_file, "%g\n", A(i)); } fclose(in_file); return 0; } int readtxtfile_vec(const char *filename, MV_Vector_double *Aptr) { FILE *in_file; char line[82]; char *line_ptr; int count=0; double tmp; MV_Vector_double &A = *Aptr; in_file = fopen( filename, "r"); if (in_file == NULL) { cerr << "Cannot open file: " << filename << endl; exit(1); } // one vector element per line while ( line_ptr = fgets(line, 82, in_file)) if (sscanf(line_ptr, "%lg", &tmp) >= 1) count++; rewind(in_file); A.newsize(count); for (int i=0; i< count; i++) { if (fscanf(in_file, "%lg", &A(i)) < 1) { printf("Error reading %s\n", filename); exit(1); } } fclose(in_file); return 0; } int writetxtfile_vec(const char *filename, const MV_Vector_int &A) { FILE *in_file; in_file = fopen(filename, "w"); if (in_file == NULL) { fprintf(stderr,"Cannot open file: %s\n", filename ); exit(1); } int N = A.size(); for (int i=0; i<N; i++) { fprintf(in_file, "%d\n", A(i)); } fclose(in_file); return 0; } int readtxtfile_vec(const char *filename, MV_Vector_int *Aptr) { FILE *in_file; char line[82]; char *line_ptr; int count=0; MV_Vector_int &A = *Aptr; int tmp; in_file = fopen( filename, "r"); if (in_file == NULL) { cerr << "Cannot open file: " << filename << endl; exit(1); } // one vector element per line while ( line_ptr = fgets(line, 82, in_file)) if (sscanf(line_ptr, "%d", &tmp) >= 1) count++; rewind(in_file); A.newsize(count); for (int i=0; i< count; i++) { if (fscanf(in_file, "%d", &A(i)) < 1) { printf("Error reading %s\n", filename); exit(1); } } fclose(in_file); return 0; }

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #include <stdio.h> #include <stdlib.h> #include "vecdefs.h" #include VECTOR_H int TestVec_double(int N, int start, int end); int TestVec_int(int N, int start, int end); int TestVec_float(int N, int start, int end); int TestBlas1_double(int N); int TestBlas1_float(int N); int TestMat_double(int M, int N, int Istart, int Jstart, int Iend, int Jend); int TestMat_int(int M, int N, int Istart, int Jstart, int Iend, int Jend); int TestMat_float(int M, int N, int Istart, int Jstart, int Iend, int Jend); int TestBlas1_COMPLEX(int N); int TestVec_COMPLEX(int N, int start, int end); int TestMat_COMPLEX(int M, int N, int Istart, int Jstart, int Iend, int Jend); int main(int argc, char *argv[]) { int M; int N; int Istart, Iend; int Jstart, Jend; if (argc < 7) { printf("Usage: M, N, Istart, Jstart, Iend, Jend\n"); exit(1); } M = atoi(argv[1]); N = atoi(argv[2]); Istart = atoi(argv[3]); Jstart = atoi(argv[4]); Iend = atoi(argv[5]); Jend = atoi(argv[6]); TestVec_double(N, Istart, Iend); TestVec_float(N, Istart, Iend); TestVec_int(N, Istart, Iend); TestMat_double(M, N, Istart, Jstart, Iend, Jend); TestMat_float(M, N, Istart, Jstart, Iend, Jend); TestMat_int(M, N, Istart, Jstart, Iend, Jend); TestBlas1_double(N); TestBlas1_float(N); TestVec_COMPLEX(N, Istart, Iend); TestMat_COMPLEX(M, N, Istart, Jstart, Iend, Jend); TestBlas1_COMPLEX(N); cout << endl << endl; cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl; cout << "+ Successful completion of testing for MV++ +" << endl; cout << "+ No errors detected in conversion or blas routines. +" << endl; cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl; cout << endl << endl; return 0; }

View file

@ -0,0 +1 @@
include ../makefile.def MV_INCLUDE_DIR = ../include INCLUDES = -I$(MV_INCLUDE_DIR) VECTOR_H = "mvvall.h" VECTOR_DOUBLE = MV_Vector_double VECTOR_FLOAT = MV_Vector_float VECTOR_INT = MV_Vector_int VECTOR_COMPLEX = MV_Vector_COMPLEX MATRIX_H = "mvm.h" MATRIX_DOUBLE = MV_ColMat_double MATRIX_FLOAT = MV_ColMat_float MATRIX_INT = MV_ColMat_int MATRIX_COMPLEX = MV_ColMat_COMPLEX # Non-templated version # SRCS = tvecd.cc \ tmatd.cc \ tblas1d.cc \ tveci.cc \ tmati.cc \ tvecf.cc \ tmatf.cc \ tblas1f.cc \ tvecc.cc \ tmatc.cc \ tblas1c.cc OBJS = tvecd.o \ tmatd.o \ tblas1d.o \ tveci.o \ tmati.o \ tvecf.o \ tmatf.o \ tblas1f.o \ tvecc.o \ tmatc.o \ tblas1c.o EXECS = main testio all: $(EXECS) main : main.o $(OBJS) $(CCC) $(CCCFLAGS) -o main main.o $(OBJS) $(LDFLAGS) testio : testio.o $(CCC) $(CCCFLAGS) -o testio testio.o $(LDFLAGS) tvecd.o : tvecd.cc \ $(MV_INCLUDE_DIR)/mvvd.h $(MV_INCLUDE_DIR)/mvvind.h tmatd.o : tmatd.cc \ $(MV_INCLUDE_DIR)/mvmd.h $(MV_INCLUDE_DIR)/mvvd.h \ $(MV_INCLUDE_DIR)/mvvind.h tblas1d.o : tblas1d.cc \ $(MV_INCLUDE_DIR)/mvmd.h $(MV_INCLUDE_DIR)/mvvd.h \ $(MV_INCLUDE_DIR)/mvvind.h tvecf.o : tvecf.cc \ $(MV_INCLUDE_DIR)/mvvf.h $(MV_INCLUDE_DIR)/mvvind.h tmatf.o : tmatf.cc \ $(MV_INCLUDE_DIR)/mvmf.h $(MV_INCLUDE_DIR)/mvvf.h \ $(MV_INCLUDE_DIR)/mvvind.h tblas1f.o : tblas1f.cc \ $(MV_INCLUDE_DIR)/mvmf.h $(MV_INCLUDE_DIR)/mvvf.h \ $(MV_INCLUDE_DIR)/mvvind.h tveci.o : tveci.cc \ $(MV_INCLUDE_DIR)/mvvi.h $(MV_INCLUDE_DIR)/mvvind.h tmati.o : tmati.cc \ $(MV_INCLUDE_DIR)/mvmi.h $(MV_INCLUDE_DIR)/mvvi.h \ $(MV_INCLUDE_DIR)/mvvind.h tvecc.o : tvecc.cc \ $(MV_INCLUDE_DIR)/mvvc.h $(MV_INCLUDE_DIR)/mvvind.h tmatc.o : tmatc.cc \ $(MV_INCLUDE_DIR)/mvmc.h $(MV_INCLUDE_DIR)/mvvc.h \ $(MV_INCLUDE_DIR)/mvvind.h tblas1c.o : tblas1c.cc \ $(MV_INCLUDE_DIR)/mvmc.h $(MV_INCLUDE_DIR)/mvvc.h \ $(MV_INCLUDE_DIR)/mvvind.h testio.o : testio.cc \ $(MV_INCLUDE_DIR)/mvmd.h $(MV_INCLUDE_DIR)/mvvi.h \ $(MV_INCLUDE_DIR)/mvvind.h wipe: clean clean: /bin/rm -f $(OBJS) $(EXECS) main.o testio.o main testio; /bin/rm -r -f ptrepository; /bin/rm -r -f Templates.DB;

View file

@ -0,0 +1,3 @@
main 8 9 2 4 6 7
testio testio.dat
rm testio0001.tmp

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // tmat.cc Testing routine for MV++ matrix class // // Usage: <M> <N> <block-index start> <block-index end> // // (Suggest these to be relatively small positive numbers, as // output is proportional to "length".) // // This test program is self-explanatory (there are diagnostics // output message of what each step is performing)and excercise // most of the MV_ColMat methods, thus providing examples of how // these are used. // // modif gŽrard Rio : // - utilisation des routines spŽcific mac pour l'entrŽe via la console // - introduction de la surcharge d'Žcriture dans la classe MVvecteur #include <iostream.h> #include <stdlib.h> #include <stdio.h> #include <console.h> // spŽcific mac : GR #include <SIOUX.h> // spŽcific mac : GR #include "mvvtp.h" int main(int argc, char *argv[]) { argc = ccommand( &argv ); // spŽcific mac : GR if (argc<2) { cout << "Usage: " << " M " << endl; exit(1); } int N = atoi(argv[1]); int j; cout << "Using N = " << N << " N = " << N << endl; MV_Vector<double> x(N), y(N), z(N); for (j=0; j<N; j++) { x(j) = j; y(j) = j / 100.0; } cout << " x " << endl; cout << x << endl; cout << "Testing MV_Vector * scalar: " << endl; z = x * 2.0; cout << " z = x * 2.0 : " << endl; for (unsigned int i=0;i< z.size();i++) cout << z(i) << " "; cout << endl; z = z * 2.0; cout << " z = z * 2.0 : " << endl; for (unsigned int i=0;i< z.size();i++) cout << z(i)<< " " ; cout << endl; cout << "Testing MV_Vector *= scalar: " << endl; (x *= 2.0) *= 2.0; cout << " (x *= 2.0) *= 2.0: " << endl; for (unsigned int i=0;i< x.size();i++) cout << x(i) << " "; cout << endl; cout << "Testing MV_Vector += MV_Vector: " << endl; x = 2.0; y = 1.0; cout << "x=2.0, y=1.0; x += y: " << endl; x += y; for (unsigned int i=0;i< x.size();i++) cout << x(i) << " "; cout << endl; cout << "Testing MV_Vector -= MV_Vector: " << endl; x = 2.0; y = 1.0; cout << "x=2.0, y=1.0; x -= y: " << endl; x -= y; for (unsigned int i=0;i< x.size();i++) cout << x(i) << " "; cout << endl; return 0; }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // tmat.cc Testing routine for MV++ matrix class // // Usage: <M> <N> <block-index start> <block-index end> // // (Suggest these to be relatively small positive numbers, as // output is proportional to "length".) // // This test program is self-explanatory (there are diagnostics // output message of what each step is performing)and excercise // most of the MV_ColMat methods, thus providing examples of how // these are used. // #include <iostream.h> #include <stdlib.h> #include "vecdefs.h" #include VECTOR_H int TestBlas1_COMPLEX(int N) { int j; cout << " \n TestBlas1_COMPLEX: \n" ; cout << "Using N = " << N << " N = " << N << endl; VECTOR_COMPLEX x(N), y(N), z(N); for (j=0; j<N; j++) { x(j) = j; y(j) = j / 100; } cout << " x " << endl; cout << x << endl; cout << "Testing Vector * scalar: " << endl; z = x * 2; cout << " z = x * 2 : " << endl; cout << z << endl; z = z * 2; cout << " z = z * 2 : " << endl; cout << z << endl; cout << "Testing Vector *= scalar: " << endl; (x *= 2) *= 2; cout << " (x *= 2) *= 2: " << endl; cout << x << endl; cout << "Testing Vector += Vector: " << endl; x = (COMPLEX) 2; y = (COMPLEX) 1; cout << "x=2, y=1; x += y: " << endl; x += y; cout << x << endl; cout << "Testing Vector -= Vector: " << endl; x = (COMPLEX) 2; y = (COMPLEX) 1; cout << "x=2, y=1; x -= y: " << endl; x -= y; cout << x << endl << endl; return 0; }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // tmat.cc Testing routine for MV++ matrix class // // Usage: <M> <N> <block-index start> <block-index end> // // (Suggest these to be relatively small positive numbers, as // output is proportional to "length".) // // This test program is self-explanatory (there are diagnostics // output message of what each step is performing)and excercise // most of the MV_ColMat methods, thus providing examples of how // these are used. // #include <iostream.h> #include <stdlib.h> #include "vecdefs.h" #include VECTOR_H int TestBlas1_double(int N) { int j; cout << " \n TestBlas1_double: \n" ; cout << "Using N = " << N << " N = " << N << endl; VECTOR_double x(N), y(N), z(N); for (j=0; j<N; j++) { x(j) = j; y(j) = j / 100; } cout << " x " << endl; cout << x << endl; cout << "Testing Vector * scalar: " << endl; z = x * 2; cout << " z = x * 2 : " << endl; cout << z << endl; z = z * 2; cout << " z = z * 2 : " << endl; cout << z << endl; cout << "Testing Vector *= scalar: " << endl; (x *= 2) *= 2; cout << " (x *= 2) *= 2: " << endl; cout << x << endl; cout << "Testing Vector += Vector: " << endl; x = (double) 2; y = (double) 1; cout << "x=2, y=1; x += y: " << endl; x += y; cout << x << endl; cout << "Testing Vector -= Vector: " << endl; x = (double) 2; y = (double) 1; cout << "x=2, y=1; x -= y: " << endl; x -= y; cout << x << endl << endl; return 0; }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // tmat.cc Testing routine for MV++ matrix class // // Usage: <M> <N> <block-index start> <block-index end> // // (Suggest these to be relatively small positive numbers, as // output is proportional to "length".) // // This test program is self-explanatory (there are diagnostics // output message of what each step is performing)and excercise // most of the MV_ColMat methods, thus providing examples of how // these are used. // #include <iostream.h> #include <stdlib.h> #include "vecdefs.h" #include VECTOR_H int TestBlas1_float(int N) { int j; cout << " \n TestBlas1_float: \n" ; cout << "Using N = " << N << " N = " << N << endl; VECTOR_float x(N), y(N), z(N); for (j=0; j<N; j++) { x(j) = j; y(j) = j / 100; } cout << " x " << endl; cout << x << endl; cout << "Testing Vector * scalar: " << endl; z = x * 2; cout << " z = x * 2 : " << endl; cout << z << endl; z = z * 2; cout << " z = z * 2 : " << endl; cout << z << endl; cout << "Testing Vector *= scalar: " << endl; (x *= 2) *= 2; cout << " (x *= 2) *= 2: " << endl; cout << x << endl; cout << "Testing Vector += Vector: " << endl; x = (float) 2; y = (float) 1; cout << "x=2, y=1; x += y: " << endl; x += y; cout << x << endl; cout << "Testing Vector -= Vector: " << endl; x = (float) 2; y = (float) 1; cout << "x=2, y=1; x -= y: " << endl; x -= y; cout << x << endl << endl; return 0; }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // tmat.cc Testing routine for MV++ matrix class // // Usage: <M> <N> <block-index start> <block-index end> // // (Suggest these to be relatively small positive numbers, as // output is proportional to "length".) // // This test program is self-explanatory (there are diagnostics // output message of what each step is performing)and excercise // most of the MV_ColMat methods, thus providing examples of how // these are used. // #include <iostream.h> #include <stdlib.h> #include "vecdefs.h" #include VECTOR_H int main(int argc, char *argv[]) { if (argc<2) { cout << "Usage: " << " M " << endl; exit(1); } int N = atoi(argv[1]); int j; cout << "Using N = " << N << " N = " << N << endl; VECTOR_int x(N), y(N), z(N); for (j=0; j<N; j++) { x(j) = j; y(j) = j / 100; } cout << " x " << endl; cout << x << endl; cout << "Testing Vector * scalar: " << endl; z = x * 2; cout << " z = x * 2 : " << endl; cout << z << endl; z = z * 2; cout << " z = z * 2 : " << endl; cout << z << endl; cout << "Testing Vector *= scalar: " << endl; (x *= 2) *= 2; cout << " (x *= 2) *= 2: " << endl; cout << x << endl; cout << "Testing Vector += Vector: " << endl; x = (int) 2; y = (int) 1; cout << "x=2, y=1; x += y: " << endl; x += y; cout << x << endl; cout << "Testing Vector -= Vector: " << endl; x = (int) 2; y = (int) 1; cout << "x=2, y=1; x -= y: " << endl; x -= y; cout << x << endl; return 0; }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // tmat.cc Testing routine for MV++ matrix class // // Usage: <M> <N> <block-index start> <block-index end> // // (Suggest these to be relatively small positive numbers, as // output is proportional to "length".) // // This test program is self-explanatory (there are diagnostics // output message of what each step is performing)and excercise // most of the MV_ColMat methods, thus providing examples of how // these are used. // #include <iostream.h> #include <stdlib.h> #include "vecdefs.h" #include VECTOR_H int TestBlas1_$TYPE(int N) { int j; cout << " \n TestBlas1_$TYPE: \n" ; cout << "Using N = " << N << " N = " << N << endl; VECTOR_$TYPE x(N), y(N), z(N); for (j=0; j<N; j++) { x(j) = j; y(j) = j / 100; } cout << " x " << endl; cout << x << endl; cout << "Testing Vector * scalar: " << endl; z = x * 2; cout << " z = x * 2 : " << endl; cout << z << endl; z = z * 2; cout << " z = z * 2 : " << endl; cout << z << endl; cout << "Testing Vector *= scalar: " << endl; (x *= 2) *= 2; cout << " (x *= 2) *= 2: " << endl; cout << x << endl; cout << "Testing Vector += Vector: " << endl; x = ($TYPE) 2; y = ($TYPE) 1; cout << "x=2, y=1; x += y: " << endl; x += y; cout << x << endl; cout << "Testing Vector -= Vector: " << endl; x = ($TYPE) 2; y = ($TYPE) 1; cout << "x=2, y=1; x -= y: " << endl; x -= y; cout << x << endl << endl; return 0; }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #include <iostream.h> #include "mvvd.h" #include "mvvi.h" #include "mvblasd.h" #include "iotext.h" const char *tmpname = "testio0001.tmp"; main(int argc, char *argv[]) { if (argc < 2) { cout << "Usage <vector.text> " << endl; exit(1); } char *filename = argv[1]; MV_Vector_double A, B; readtxtfile_vec(filename, &A); writetxtfile_vec(tmpname, A); readtxtfile_vec(filename, &B); if (norm(A-B) / norm(A+B) < 0.0000001 ) cout << "testio passed. " << endl; else cout << "testio failed.\n" << endl; }

File diff suppressed because one or more lines are too long

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ #include <iostream.h> #include <stdlib.h> #include <console.h> // spŽcific mac : GR #include <SIOUX.h> // spŽcific mac : GR #include "mvmtp.h" void foo(MV_ColMat<double> *Aptr) { MV_ColMat<double> &A = *Aptr; A = 0.0; } int main(int argc, char *argv[]) { argc = ccommand( &argv ); // spŽcific mac : GR if (argc<7) { cout << "Usage " << argv[0] << " M N start0 end0 start1 end1" << endl; cout << "\n plus clair : " << "\n la taille en ligne et colonne (peuvent <20>tre diffŽrentes -> 2 nombres " << "\n deux index de vecteur (plage d'indice) -> 4 nombres (min-max min-max" << endl; exit(1); } int M = atoi(argv[1]); int N = atoi(argv[2]); int s0 = atoi(argv[3]); int e0 = atoi(argv[4]); int s1 = atoi(argv[5]); int e1 = atoi(argv[6]); cout << "Using M = " << M << endl; cout << "Using N = " << N << endl; cout << "Using start0 = " << s0 << endl; cout << "Using end0 = " << e0 << endl; cout << "Using start1 = " << s1 << endl; cout << "Using end1 = " << e1 << endl << endl;; MV_ColMat<double> A(M,N); MV_ColMat<double> B(M,N); A = 3.0; B = 7.0; cout << "\n crŽation de A = 3 : \n" << A << " " << "\n puis de B = 7 : \n" << B << endl; A( M/2, N/2) = 0.0; cout << " A=3.0; A(M/2, N/2) = 0.0 " << endl; cout << A << endl; MV_VecIndex I(s0, e0); MV_VecIndex J(s1, e1); cout << "I = (" << s0 << ":"<< e0 <<") " << endl; cout << "J = (" << s1 << ":"<< e1 <<") " << endl; A(I,J+1) = B(I,J); cout << "A(I,J+1) = B(I,J) " << endl << A << endl; A(I,J) = 8.0; cout << "A(I,J) = 8.0 " << endl << A << endl; #if 0 // The following test section attempts to take the address of a temporary // result and hence not compilable on all C++ platforms. cout << "calling foo(&A(I)) " << endl; foo(&A(I)); cout << "A(I)" << endl << A << endl; MV_ColMat<int> C(M, N); MV_ColMat<int> &c = C(I,J); C = 1; c = 0; cout << "After renaming c to be C(I,J), and executing c = 0.0 " << endl << C << endl; #endif return (0); }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // tmat.cc Testing routine for MV++ matrix class // // Usage: <M> <N> <block-index start> <block-index end> // // (Suggest these to be relatively small positive numbers, as // output is proportional to "length".) // // This test program is self-explanatory (there are diagnostics // output message of what each step is performing)and excercise // most of the MV_ColMat methods, thus providing examples of how // these are used. // #include <iostream.h> #include <stdlib.h> #include "vecdefs.h" #include MATRIX_H // #include "systime.c" // example function to zero-out a matrix: illustrates how to pass // matrixes by "reference" (modifiable), as opposed to by "value" // void foo(MATRIX_COMPLEX *Aptr) { MATRIX_COMPLEX &A = *Aptr; A = (COMPLEX) 0; // equivalent to "Aptr->operator=(0);" } // example funciton to return MV_ColMats: (returns matrix + scalar) // MATRIX_COMPLEX fexample(const MATRIX_COMPLEX A, COMPLEX s) { int M = A.size(0); int N = A.size(1); MATRIX_COMPLEX B(M,N); for (int j=0; j<N; j++) for (int i=0; i<M; i++) B(i,j) = A(i,j) + s; return B; } int TestMat_COMPLEX(int M, int N, int Istart, int Jstart, int Iend, int Jend) { int i,j; cout << endl; cout << "TestMat_COMPLEX : " << endl; cout << "Using M = " << M << " N = " << N << " Istart = " << Istart << " Iend = " << Iend << " Jstart = " << Jstart << " Jend = " << Jend << endl; MATRIX_COMPLEX C; cout << " Test null constructor: MV_ColMat C() " << endl; cout << C << endl; MATRIX_COMPLEX A(M,N); cout << " Test MV_ColMat(int,int) constructor: MV_ColMat A(M,N) " << endl; cout << " values should be uninitalized.. " << endl; cout << A << endl; MATRIX_COMPLEX B(M,N, 3); cout << " Test MV_ColMat(int, int, val) constructor: MV_ColMat B(N, 3.0) " << endl; cout << " all values should be 3: " << endl; cout << B << endl; cout << " Test A(i,j) indexing, set A(i,j) = 100*j + i, and B(i) = -A(i,j)" << endl; for (j=0; j<N; j++) for (i=0; i<M; i++) { A(i,j) = 100*i+j; B(i,j) = - A(i,j); } cout << " A " << endl; cout << A << endl; cout << " B " << endl; cout << B << endl; cout << "Testing MV_VecIndex I(Istart, Iend) " << endl; // MV_VecIndex I(Istart, Iend); MV_VecIndex J(Jstart, Jend); // cout << "I = (" << Istart << ":"<< Iend <<") " << endl; cout << "J = (" << Jstart << ":"<< Jend <<") " << endl; cout << "Test A(I,J) = B(I,J) " << endl ; A(I,J) = B(I,J); cout << A << endl; cout << "Testing A(I,J) = 1.1 " << endl; A(I,J) = 11; cout << "A" << endl << A << endl; cout << "Testing MV_VecIndex + operators: " << endl; cout << " A = 0.0; A(I+1,J-1) = 1.0; " << endl; // A = (COMPLEX) 0; A(I+1,J-1) = (COMPLEX) 1; // cout << " A " << endl << A << endl; cout << endl; return 0; }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // tmat.cc Testing routine for MV++ matrix class // // Usage: <M> <N> <block-index start> <block-index end> // // (Suggest these to be relatively small positive numbers, as // output is proportional to "length".) // // This test program is self-explanatory (there are diagnostics // output message of what each step is performing)and excercise // most of the MV_ColMat methods, thus providing examples of how // these are used. // #include <iostream.h> #include <stdlib.h> #include "vecdefs.h" #include MATRIX_H // #include "systime.c" // example function to zero-out a matrix: illustrates how to pass // matrixes by "reference" (modifiable), as opposed to by "value" // void foo(MATRIX_double *Aptr) { MATRIX_double &A = *Aptr; A = (double) 0; // equivalent to "Aptr->operator=(0);" } // example funciton to return MV_ColMats: (returns matrix + scalar) // MATRIX_double fexample(const MATRIX_double A, double s) { int M = A.size(0); int N = A.size(1); MATRIX_double B(M,N); for (int j=0; j<N; j++) for (int i=0; i<M; i++) B(i,j) = A(i,j) + s; return B; } int TestMat_double(int M, int N, int Istart, int Jstart, int Iend, int Jend) { int i,j; cout << endl; cout << "TestMat_double : " << endl; cout << "Using M = " << M << " N = " << N << " Istart = " << Istart << " Iend = " << Iend << " Jstart = " << Jstart << " Jend = " << Jend << endl; MATRIX_double C; cout << " Test null constructor: MV_ColMat C() " << endl; cout << C << endl; MATRIX_double A(M,N); cout << " Test MV_ColMat(int,int) constructor: MV_ColMat A(M,N) " << endl; cout << " values should be uninitalized.. " << endl; cout << A << endl; MATRIX_double B(M,N, 3); cout << " Test MV_ColMat(int, int, val) constructor: MV_ColMat B(N, 3.0) " << endl; cout << " all values should be 3: " << endl; cout << B << endl; cout << " Test A(i,j) indexing, set A(i,j) = 100*j + i, and B(i) = -A(i,j)" << endl; for (j=0; j<N; j++) for (i=0; i<M; i++) { A(i,j) = 100*i+j; B(i,j) = - A(i,j); } cout << " A " << endl; cout << A << endl; cout << " B " << endl; cout << B << endl; cout << "Testing MV_VecIndex I(Istart, Iend) " << endl; // MV_VecIndex I(Istart, Iend); MV_VecIndex J(Jstart, Jend); // cout << "I = (" << Istart << ":"<< Iend <<") " << endl; cout << "J = (" << Jstart << ":"<< Jend <<") " << endl; cout << "Test A(I,J) = B(I,J) " << endl ; A(I,J) = B(I,J); cout << A << endl; cout << "Testing A(I,J) = 1.1 " << endl; A(I,J) = 11; cout << "A" << endl << A << endl; cout << "Testing MV_VecIndex + operators: " << endl; cout << " A = 0.0; A(I+1,J-1) = 1.0; " << endl; // A = (double) 0; A(I+1,J-1) = (double) 1; // cout << " A " << endl << A << endl; cout << endl; return 0; }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // tmat.cc Testing routine for MV++ matrix class // // Usage: <M> <N> <block-index start> <block-index end> // // (Suggest these to be relatively small positive numbers, as // output is proportional to "length".) // // This test program is self-explanatory (there are diagnostics // output message of what each step is performing)and excercise // most of the MV_ColMat methods, thus providing examples of how // these are used. // #include <iostream.h> #include <stdlib.h> #include "vecdefs.h" #include MATRIX_H // #include "systime.c" // example function to zero-out a matrix: illustrates how to pass // matrixes by "reference" (modifiable), as opposed to by "value" // void foo(MATRIX_float *Aptr) { MATRIX_float &A = *Aptr; A = (float) 0; // equivalent to "Aptr->operator=(0);" } // example funciton to return MV_ColMats: (returns matrix + scalar) // MATRIX_float fexample(const MATRIX_float A, float s) { int M = A.size(0); int N = A.size(1); MATRIX_float B(M,N); for (int j=0; j<N; j++) for (int i=0; i<M; i++) B(i,j) = A(i,j) + s; return B; } int TestMat_float(int M, int N, int Istart, int Jstart, int Iend, int Jend) { int i,j; cout << endl; cout << "TestMat_float : " << endl; cout << "Using M = " << M << " N = " << N << " Istart = " << Istart << " Iend = " << Iend << " Jstart = " << Jstart << " Jend = " << Jend << endl; MATRIX_float C; cout << " Test null constructor: MV_ColMat C() " << endl; cout << C << endl; MATRIX_float A(M,N); cout << " Test MV_ColMat(int,int) constructor: MV_ColMat A(M,N) " << endl; cout << " values should be uninitalized.. " << endl; cout << A << endl; MATRIX_float B(M,N, 3); cout << " Test MV_ColMat(int, int, val) constructor: MV_ColMat B(N, 3.0) " << endl; cout << " all values should be 3: " << endl; cout << B << endl; cout << " Test A(i,j) indexing, set A(i,j) = 100*j + i, and B(i) = -A(i,j)" << endl; for (j=0; j<N; j++) for (i=0; i<M; i++) { A(i,j) = 100*i+j; B(i,j) = - A(i,j); } cout << " A " << endl; cout << A << endl; cout << " B " << endl; cout << B << endl; cout << "Testing MV_VecIndex I(Istart, Iend) " << endl; // MV_VecIndex I(Istart, Iend); MV_VecIndex J(Jstart, Jend); // cout << "I = (" << Istart << ":"<< Iend <<") " << endl; cout << "J = (" << Jstart << ":"<< Jend <<") " << endl; cout << "Test A(I,J) = B(I,J) " << endl ; A(I,J) = B(I,J); cout << A << endl; cout << "Testing A(I,J) = 1.1 " << endl; A(I,J) = 11; cout << "A" << endl << A << endl; cout << "Testing MV_VecIndex + operators: " << endl; cout << " A = 0.0; A(I+1,J-1) = 1.0; " << endl; // A = (float) 0; A(I+1,J-1) = (float) 1; // cout << " A " << endl << A << endl; cout << endl; return 0; }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // tmat.cc Testing routine for MV++ matrix class // // Usage: <M> <N> <block-index start> <block-index end> // // (Suggest these to be relatively small positive numbers, as // output is proportional to "length".) // // This test program is self-explanatory (there are diagnostics // output message of what each step is performing)and excercise // most of the MV_ColMat methods, thus providing examples of how // these are used. // #include <iostream.h> #include <stdlib.h> #include "vecdefs.h" #include MATRIX_H // #include "systime.c" // example function to zero-out a matrix: illustrates how to pass // matrixes by "reference" (modifiable), as opposed to by "value" // void foo(MATRIX_int *Aptr) { MATRIX_int &A = *Aptr; A = (int) 0; // equivalent to "Aptr->operator=(0);" } // example funciton to return MV_ColMats: (returns matrix + scalar) // MATRIX_int fexample(const MATRIX_int A, int s) { int M = A.size(0); int N = A.size(1); MATRIX_int B(M,N); for (int j=0; j<N; j++) for (int i=0; i<M; i++) B(i,j) = A(i,j) + s; return B; } int TestMat_int(int M, int N, int Istart, int Jstart, int Iend, int Jend) { int i,j; cout << endl; cout << "TestMat_int : " << endl; cout << "Using M = " << M << " N = " << N << " Istart = " << Istart << " Iend = " << Iend << " Jstart = " << Jstart << " Jend = " << Jend << endl; MATRIX_int C; cout << " Test null constructor: MV_ColMat C() " << endl; cout << C << endl; MATRIX_int A(M,N); cout << " Test MV_ColMat(int,int) constructor: MV_ColMat A(M,N) " << endl; cout << " values should be uninitalized.. " << endl; cout << A << endl; MATRIX_int B(M,N, 3); cout << " Test MV_ColMat(int, int, val) constructor: MV_ColMat B(N, 3.0) " << endl; cout << " all values should be 3: " << endl; cout << B << endl; cout << " Test A(i,j) indexing, set A(i,j) = 100*j + i, and B(i) = -A(i,j)" << endl; for (j=0; j<N; j++) for (i=0; i<M; i++) { A(i,j) = 100*i+j; B(i,j) = - A(i,j); } cout << " A " << endl; cout << A << endl; cout << " B " << endl; cout << B << endl; cout << "Testing MV_VecIndex I(Istart, Iend) " << endl; // MV_VecIndex I(Istart, Iend); MV_VecIndex J(Jstart, Jend); // cout << "I = (" << Istart << ":"<< Iend <<") " << endl; cout << "J = (" << Jstart << ":"<< Jend <<") " << endl; cout << "Test A(I,J) = B(I,J) " << endl ; A(I,J) = B(I,J); cout << A << endl; cout << "Testing A(I,J) = 1.1 " << endl; A(I,J) = 11; cout << "A" << endl << A << endl; cout << "Testing MV_VecIndex + operators: " << endl; cout << " A = 0.0; A(I+1,J-1) = 1.0; " << endl; // A = (int) 0; A(I+1,J-1) = (int) 1; // cout << " A " << endl << A << endl; cout << endl; return 0; }

View file

@ -0,0 +1 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /* */ /* */ /* MV++ Numerical Matrix/Vector C++ Library */ /* MV++ Version 1.5 */ /* */ /* R. Pozo */ /* National Institute of Standards and Technology */ /* */ /* NOTICE */ /* */ /* Permission to use, copy, modify, and distribute this software and */ /* its documentation for any purpose and without fee is hereby granted */ /* provided that this permission notice appear in all copies and */ /* supporting documentation. */ /* */ /* Neither the Institution (National Institute of Standards and Technology) */ /* nor the author makes any representations about the suitability of this */ /* software for any purpose. This software is provided ``as is''without */ /* expressed or implied warranty. */ /* */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ // // tmat.cc Testing routine for MV++ matrix class // // Usage: <M> <N> <block-index start> <block-index end> // // (Suggest these to be relatively small positive numbers, as // output is proportional to "length".) // // This test program is self-explanatory (there are diagnostics // output message of what each step is performing)and excercise // most of the MV_ColMat methods, thus providing examples of how // these are used. // #include <iostream.h> #include <stdlib.h> #include "vecdefs.h" #include MATRIX_H // #include "systime.c" // example function to zero-out a matrix: illustrates how to pass // matrixes by "reference" (modifiable), as opposed to by "value" // void foo(MATRIX_$TYPE *Aptr) { MATRIX_$TYPE &A = *Aptr; A = ($TYPE) 0; // equivalent to "Aptr->operator=(0);" } // example funciton to return MV_ColMats: (returns matrix + scalar) // MATRIX_$TYPE fexample(const MATRIX_$TYPE A, $TYPE s) { int M = A.size(0); int N = A.size(1); MATRIX_$TYPE B(M,N); for (int j=0; j<N; j++) for (int i=0; i<M; i++) B(i,j) = A(i,j) + s; return B; } int TestMat_$TYPE(int M, int N, int Istart, int Jstart, int Iend, int Jend) { int i,j; cout << endl; cout << "TestMat_$TYPE : " << endl; cout << "Using M = " << M << " N = " << N << " Istart = " << Istart << " Iend = " << Iend << " Jstart = " << Jstart << " Jend = " << Jend << endl; MATRIX_$TYPE C; cout << " Test null constructor: MV_ColMat C() " << endl; cout << C << endl; MATRIX_$TYPE A(M,N); cout << " Test MV_ColMat(int,int) constructor: MV_ColMat A(M,N) " << endl; cout << " values should be uninitalized.. " << endl; cout << A << endl; MATRIX_$TYPE B(M,N, 3); cout << " Test MV_ColMat(int, int, val) constructor: MV_ColMat B(N, 3.0) " << endl; cout << " all values should be 3: " << endl; cout << B << endl; cout << " Test A(i,j) indexing, set A(i,j) = 100*j + i, and B(i) = -A(i,j)" << endl; for (j=0; j<N; j++) for (i=0; i<M; i++) { A(i,j) = 100*i+j; B(i,j) = - A(i,j); } cout << " A " << endl; cout << A << endl; cout << " B " << endl; cout << B << endl; cout << "Testing MV_VecIndex I(Istart, Iend) " << endl; // MV_VecIndex I(Istart, Iend); MV_VecIndex J(Jstart, Jend); // cout << "I = (" << Istart << ":"<< Iend <<") " << endl; cout << "J = (" << Jstart << ":"<< Jend <<") " << endl; cout << "Test A(I,J) = B(I,J) " << endl ; A(I,J) = B(I,J); cout << A << endl; cout << "Testing A(I,J) = 1.1 " << endl; A(I,J) = 11; cout << "A" << endl << A << endl; cout << "Testing MV_VecIndex + operators: " << endl; cout << " A = 0.0; A(I+1,J-1) = 1.0; " << endl; // A = ($TYPE) 0; A(I+1,J-1) = ($TYPE) 1; // cout << " A " << endl << A << endl; cout << endl; return 0; }

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1 @@
SparseLib++ v. 1.5c: Numerical Sparse Matrix Classes in C++ SparseLib++ is a C++ class library for efficient sparse matrix computations across various computational platforms. The software package consists of matrix objects representing several sparse storage formats currently in use (in this release: compressed row, compressed column and coordinate formats), providing basic functionality for managing sparse matrices, together with efficient kernel mathematical operations (e.g. sparse matrix-vector multiply). Routines based on the Sparse BLAS are used to enhance portability and performance. Included in the package are various preconditioners commonly used in iterative solvers for linear systems of equations. The focus is on computational support for iterative methods, but the sparse matrix objects presented here can be used on their own. 0) What's new in v. 1.5 ----------------------- o) supports ANSI C++ templated complex classes (e.g. complex<double>) for GNU g++ 2.7.0 o) uses new MV++ vector/matrix 1.3 classes o) several bug fixes for reading Harwell-Boeing files and converting between matrix storage types o) supports 'P' format specifier in Harwell-Boeing files o) Bug fixes in v. 1.5c i) Added check for end-of row/column in elimination loop of ilupre.cc and ilupre_double.cc ii) fixed "nz_ = nz_" bug in coord_double::newsize() 1) More About SparseLib++ ------------------------ i) User's Guide and man-style pages are available via WWW: http://gams.cam.nist.gov/acmd/Staff/RPozo/sparselib++.html ftp: gams.nist.gov:~ftp/pub/pozo/docs/sparselib++.ps.Z ii) Code examples are in ./testing/*.cc 2) Installing SparseLib++ library ---------------------------------- You'll need to build the following three libraries ./lib/libmv.a MV++ basic matrix/vector library ./lib/libsparselib.a sparse matrix library ./lib/libspblas.a sparse BLAS library i) cd to root directory where SparseLib++ will be installed ii) edit makefile.def to specify your specify your particular C++ compiler iii) type "make sp". ("make" by itself will provide a list of options.) 3) Testing SparseLib++ ---------------------- i) cd to SparseLib++ root dir ii) "make test" will run a test suites and leave their output in ./testing/sp_test.out. 4) Package components --------------------- i) SparseLib++ and Sparse BLAS include files are in ./include ii) SparseLib++ and Sparse BLAS libraries are in ./lib iii) Lower-level Fortran-compatible Matrix/Vector library (MV++)is in ./lib/libmv.a 5) Help! -------- Questions, comments, suggestions, etc. can be sent to pozo@cam.nist.gov. 6) #include <std/disclaimer.h> ------------------------------ Be aware that SparseLib++ has been tested with GNU g++ 2.6.3, and 2.7.0 together with Sun C++ (CC) v. 4.0.1. and may not necessarily work with earlier versions of these compilers. To find out the version of your compiler use "g++ -v" or "CC -V foo".

View file

@ -0,0 +1,126 @@
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* ******** *** SparseLib++ */
/* ******* ** *** *** *** v. 1.5c */
/* ***** *** ******** ******** */
/* ***** *** ******** ******** R. Pozo */
/* ** ******* *** ** *** *** K. Remington */
/* ******** ******** A. Lumsdaine */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* */
/* */
/* SparseLib++ : Sparse Matrix Library */
/* */
/* National Institute of Standards and Technology */
/* University of Notre Dame */
/* Authors: R. Pozo, K. Remington, A. Lumsdaine */
/* */
/* NOTICE */
/* */
/* Permission to use, copy, modify, and distribute this software and */
/* its documentation for any purpose and without fee is hereby granted */
/* provided that the above notice appear in all copies and supporting */
/* documentation. */
/* */
/* Neither the Institutions (National Institute of Standards and Technology, */
/* University of Notre Dame) nor the Authors make any representations about */
/* the suitability of this software for any purpose. This software is */
/* provided ``as is'' without expressed or implied warranty. */
/* */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* Compressed column sparse matrix (O-based, Fortran) */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
#ifndef CompCol_Mat_double_H
#define CompCol_Mat_double_H
//#include "vecdefs.h"
#include "vecdefs_GR.h" // modif GR
//#include VECTOR_H // modif GR
class CompRow_Mat_double;
class Coord_Mat_double;
class CompCol_Mat_double {
// private:
protected : // modif GR
VECTOR_double val_; // data values (nz_ elements)
VECTOR_int rowind_; // row_ind (nz_ elements)
VECTOR_int colptr_; // col_ptr (dim_[1]+1 elements)
int base_; // index base: offset of first element
int nz_; // number of nonzeros
int dim_[2]; // number of rows, cols
public:
CompCol_Mat_double(void);
CompCol_Mat_double(const CompCol_Mat_double &S);
CompCol_Mat_double(const CompRow_Mat_double &R);
CompCol_Mat_double(const Coord_Mat_double &CO);
// constructeur parmettant la création d'une matrice creuse
// M : nombre de lignes
// N : nombre de colonnes
// nz : le nombre de valeurs non nulles
// val : le vecteur qui contiend les données non nulles, sous forme
// de la suite des colonnes
// r : la suite des numéros de position des valeurs non nulles
// dans les différentes colonnes
// c : les numéros dans r de chaque début de colonne + le numéro
// du dernier élément c-a-d nz
// base : indique si la matrice est une copie ou non d'une autre matrice
CompCol_Mat_double(int M, int N, int nz, double *val, int *r,
int *c, int base=0);
CompCol_Mat_double(int M, int N, int nz, const VECTOR_double &val,
const VECTOR_int &r, const VECTOR_int &c,
int base=0);
~CompCol_Mat_double() {};
/*******************************/
/* Access and info functions */
/*******************************/
double& val(int i) { return val_(i); }
int& row_ind(int i) { return rowind_(i); }
int& col_ptr(int i) { return colptr_(i);}
const double& val(int i) const { return val_(i); }
const int& row_ind(int i) const { return rowind_(i); }
const int& col_ptr( int i) const { return (colptr_(i));}
int dim(int i) const {return dim_[i];};
int size(int i) const {return dim_[i];};
int NumNonzeros() const {return nz_;};
int base() const {return base_;}
/*******************************/
/* Assignment operator */
/*******************************/
CompCol_Mat_double& operator=(const CompCol_Mat_double &C);
CompCol_Mat_double& newsize(int M, int N, int nz);
/***********************************/
/* General access function (slow) */
/***********************************/
double operator() (int i, int j) const;
double& set(int i, int j);
/* ajout GR */
bool Exista(int i, int j) const;
/***********************************/
/* Matrix/Vector multiply */
/***********************************/
VECTOR_double operator*(const VECTOR_double &x) const;
VECTOR_double trans_mult(const VECTOR_double &x) const;
};
ostream& operator << (ostream & os, const CompCol_Mat_double & mat);
void readHB(const char *, CompCol_Mat_double &);
#endif
/* CompCol_Mat_double_H */

File diff suppressed because one or more lines are too long

Some files were not shown because too many files have changed in this diff Show more