diff --git a/src/libintegratorSparse/centralDifferencesSparse.h b/src/libintegratorSparse/centralDifferencesSparse.h
index f60da01710289cd1dd01f24577a6e1fb3b14ea4d..1b616224ba9287032a3e0dcb229bb56106547637 100644
--- a/src/libintegratorSparse/centralDifferencesSparse.h
+++ b/src/libintegratorSparse/centralDifferencesSparse.h
@@ -28,7 +28,7 @@
/*
-A class to timestep dynamics (e.g. nonlinear deformable FEM)
+A class to timestep dynamics (e.g. nonlinear deformable FEM)
using the EXPLICIT central differences integrator. Normally, the implicit newmark and
implicit backward Euler classes are recommended for nonlinear deformable FEM
due to stability under large deformations. See also integratorBase.h.
@@ -37,7 +37,7 @@ This implementation follows
WRIGGERS P.: Computational Contact Mechanics. John
Wiley & Sons, Ltd., 2002., page 275
-This class supports two damping modes:
+This class supports two damping modes:
1. Standard Rayleigh model
D = dampingMassCoef * M + dampingStiffnessCoef * K
where K is the stiffness matrix AT THE ORIGIN
@@ -45,16 +45,16 @@ This class supports two damping modes:
D = dampingMassCoef * M + dampingStiffnessCoef * K(q)
where K is the tangential stiffness matrix in the CURRENT deformed configuration q
-Mode 1. is computationally faster as it does not require system updates (i.e.,
-matrix inversion). Mode 1. is useful, for example, if you want to timestep
+Mode 1. is computationally faster as it does not require system updates (i.e.,
+matrix inversion). Mode 1. is useful, for example, if you want to timestep
linear modal analysis simulations (stiffness matrix is constant in that case).
Mode 2. gives a better damping model for large deformations, but because the
system matrix changes, requires factoring a linear system anew at each timestep.
-In order to use this class, you need to set the timestep very small, or
-else the explicit integrator will go unstable. Roughly speaking, the timestep
-must resolve the highest frequency present in your simulation.
+In order to use this class, you need to set the timestep very small, or
+else the explicit integrator will go unstable. Roughly speaking, the timestep
+must resolve the highest frequency present in your simulation.
See also integratorBase.h .
@@ -86,13 +86,13 @@ public:
inline virtual void SetTimestep(double timestep) { this->timestep = timestep; DecomposeSystemMatrix(); }
// performs one timestep of simulation
- virtual int DoTimestep();
+ virtual int DoTimestep();
- // sets q, and (optionally) qvel
- // returns 0
+ // sets q, and (optionally) qvel
+ // returns 0
virtual int SetState(double * q, double * qvel=NULL);
- // tangentialDampingMode:
+ // tangentialDampingMode:
// 0 = no updates of the damping matrix under deformations (not recommended for large deformations)
// 1 = update at every timestep (default)
// k>1 = update every kth timestep
diff --git a/src/libintegratorSparse/eulerSparse.cpp b/src/libintegratorSparse/eulerSparse.cpp
index 6d4b12e472a4d8ff6b892cc307ab3485c83d6dcd..7e0d8d0e9a6053bbc9039c108b28d0903d903fa0 100644
--- a/src/libintegratorSparse/eulerSparse.cpp
+++ b/src/libintegratorSparse/eulerSparse.cpp
@@ -97,7 +97,7 @@ int EulerSparse::DoTimestep()
// store current state
for(int i=0; i<r; i++)
{
- q_1[i] = q[i];
+ q_1[i] = q[i];
qvel_1[i] = qvel[i];
qaccel_1[i] = qaccel[i];
}
@@ -167,7 +167,7 @@ int EulerSparse::DoTimestep()
{
qvel[i] += timestep * qdelta[i];
q[i] += timestep * qvel[i];
- }
+ }
}
else
{
diff --git a/src/libintegratorSparse/eulerSparse.h b/src/libintegratorSparse/eulerSparse.h
index 31ae91dea93160a88599e5a1ddc1ea366ca4ae32..1021090e5513bdfb1102cc409d13b8ef57b0095f 100644
--- a/src/libintegratorSparse/eulerSparse.h
+++ b/src/libintegratorSparse/eulerSparse.h
@@ -66,15 +66,15 @@ public:
virtual ~EulerSparse();
- // sets q, and (optionally) qvel
- // returns 0
+ // sets q, and (optionally) qvel
+ // returns 0
virtual int SetState(double * q, double * qvel=NULL);
- virtual int DoTimestep();
+ virtual int DoTimestep();
protected:
int symplectic;
-
+
#ifdef PARDISO
PardisoSolver * pardisoSolver;
#endif
diff --git a/src/libintegratorSparse/implicitBackwardEulerSparse.cpp b/src/libintegratorSparse/implicitBackwardEulerSparse.cpp
index cf0fbf12c9220e1c49f94a0b193e75b8b706678b..84281053ec5c74736517d25e371c159b74533ee7 100644
--- a/src/libintegratorSparse/implicitBackwardEulerSparse.cpp
+++ b/src/libintegratorSparse/implicitBackwardEulerSparse.cpp
@@ -70,7 +70,7 @@ int ImplicitBackwardEulerSparse::DoTimestep()
for(int i=0; i<r; i++)
{
qaccel_1[i] = qaccel[i] = 0; // acceleration is actually not used in this integrator
- q_1[i] = q[i];
+ q_1[i] = q[i];
qvel_1[i] = qvel[i];
}
@@ -121,7 +121,7 @@ int ImplicitBackwardEulerSparse::DoTimestep()
tangentStiffnessMatrix->ScalarMultiply(dampingStiffnessCoef, rayleighDampingMatrix);
rayleighDampingMatrix->AddSubMatrix(dampingMassCoef, *massMatrix);
- // build effective stiffness:
+ // build effective stiffness:
// Keff = M + h D + h^2 * K
// compute force residual, store it into aux variable qresidual
// qresidual = h * (-D qdot - fint + fext - h * K * qdot)) // this is semi-implicit Euler
@@ -143,7 +143,7 @@ int ImplicitBackwardEulerSparse::DoTimestep()
tangentStiffnessMatrix->MultiplyVectorAdd(qvel, qresidual);
*tangentStiffnessMatrix *= timestep;
tangentStiffnessMatrix->AddSubMatrix(1.0, *massMatrix);
-
+
// add externalForces, internalForces
for(int i=0; i<r; i++)
{
@@ -192,7 +192,7 @@ int ImplicitBackwardEulerSparse::DoTimestep()
//printf("numIter: %d error2: %G\n", numIter, error);
// on the first iteration, compute initial error
- if (numIter == 0)
+ if (numIter == 0)
{
error0 = error;
errorQuotient = 1.0;
@@ -200,7 +200,7 @@ int ImplicitBackwardEulerSparse::DoTimestep()
else
{
// error divided by the initial error, before performing this iteration
- errorQuotient = error / error0;
+ errorQuotient = error / error0;
}
if (errorQuotient < epsilon * epsilon)
diff --git a/src/libintegratorSparse/implicitBackwardEulerSparse.h b/src/libintegratorSparse/implicitBackwardEulerSparse.h
index f6ead0f8c86cff8553a9e3ec520983f2a51a4ccf..3c21088ebd101e2d66d6e26bfe2380863e92a2ce 100644
--- a/src/libintegratorSparse/implicitBackwardEulerSparse.h
+++ b/src/libintegratorSparse/implicitBackwardEulerSparse.h
@@ -42,14 +42,14 @@ public:
// constrainedDOFs is an integer array of degrees of freedom that are to be fixed to zero (e.g., to permanently fix a vertex in a deformable simulation)
// constrainedDOFs are 0-indexed (separate DOFs for x,y,z), and must be pre-sorted (ascending)
// numThreads applies only to the PARDISO and SPOOLES solvers; if numThreads > 0, the sparse linear solves are multi-threaded; default: 0 (use single-threading)
- ImplicitBackwardEulerSparse(int r, double timestep, SparseMatrix * massMatrix, ForceModel * forceModel, int positiveDefiniteSolver=0, int numConstrainedDOFs=0, int * constrainedDOFs=NULL, double dampingMassCoef=0.0, double dampingStiffnessCoef=0.0, int maxIterations = 1, double epsilon = 1E-6, int numSolverThreads=0);
+ ImplicitBackwardEulerSparse(int r, double timestep, SparseMatrix * massMatrix, ForceModel * forceModel, int positiveDefiniteSolver=0, int numConstrainedDOFs=0, int * constrainedDOFs=NULL, double dampingMassCoef=0.0, double dampingStiffnessCoef=0.0, int maxIterations = 1, double epsilon = 1E-6, int numSolverThreads=0);
virtual ~ImplicitBackwardEulerSparse();
- // sets q, and (optionally) qvel
- // returns 0
+ // sets q, and (optionally) qvel
+ // returns 0
virtual int SetState(double * q, double * qvel=NULL);
- virtual int DoTimestep();
+ virtual int DoTimestep();
protected:
};
diff --git a/src/libintegratorSparse/implicitNewmarkSparse.cpp b/src/libintegratorSparse/implicitNewmarkSparse.cpp
index 7bf583c9b3d0f16c14cac1e8c66cb9f06ef3aa70..6f1b3dd1a53bad676eca00725821f631fc80df78 100644
--- a/src/libintegratorSparse/implicitNewmarkSparse.cpp
+++ b/src/libintegratorSparse/implicitNewmarkSparse.cpp
@@ -37,7 +37,7 @@
ImplicitNewmarkSparse::ImplicitNewmarkSparse(int r, double timestep, SparseMatrix * massMatrix_, ForceModel * forceModel_, int positiveDefiniteSolver_, int numConstrainedDOFs_, int * constrainedDOFs_, double dampingMassCoef, double dampingStiffnessCoef, int maxIterations, double epsilon, double NewmarkBeta, double NewmarkGamma, int numSolverThreads_): IntegratorBaseSparse(r, timestep, massMatrix_, forceModel_, numConstrainedDOFs_, constrainedDOFs_, dampingMassCoef, dampingStiffnessCoef), positiveDefiniteSolver(positiveDefiniteSolver_), numSolverThreads(numSolverThreads_)
{
this->maxIterations = maxIterations; // maxIterations = 1 for semi-implicit
- this->epsilon = epsilon;
+ this->epsilon = epsilon;
this->NewmarkBeta = NewmarkBeta;
this->NewmarkGamma = NewmarkGamma;
@@ -117,7 +117,7 @@ int ImplicitNewmarkSparse::SetState(double * q_, double * qvel_)
for(int i=0; i<numConstrainedDOFs; i++)
q[constrainedDOFs[i]] = qvel[constrainedDOFs[i]] = 0.0;
- // M * qaccel + C * qvel + R(q) = P_0
+ // M * qaccel + C * qvel + R(q) = P_0
// R(q) = P_0 = 0
// i.e. M * qaccel = - C * qvel - R(q)
@@ -171,12 +171,12 @@ int ImplicitNewmarkSparse::SetState(double * q_, double * qvel_)
printf("Error: %s sparse solver returned non-zero exit status %d.\n", solverString, (int)info);
return 1;
}
-
+
InsertRows(r, buffer, qaccel, numConstrainedDOFs, constrainedDOFs);
return 0;
}
-
+
int ImplicitNewmarkSparse::DoTimestep()
{
int numIter = 0;
@@ -187,7 +187,7 @@ int ImplicitNewmarkSparse::DoTimestep()
// store current amplitudes and set initial guesses for qaccel, qvel
for(int i=0; i<r; i++)
{
- q_1[i] = q[i];
+ q_1[i] = q[i];
qvel_1[i] = qvel[i];
qaccel_1[i] = qaccel[i];
@@ -242,7 +242,7 @@ int ImplicitNewmarkSparse::DoTimestep()
tangentStiffnessMatrix->AddSubMatrix(alpha4, *dampingMatrix, 1);
tangentStiffnessMatrix->AddSubMatrix(alpha1, *massMatrix);
-
+
// compute force residual, store it into aux variable qresidual
// qresidual = M * qaccel + C * qvel - externalForces + internalForces
@@ -281,7 +281,7 @@ int ImplicitNewmarkSparse::DoTimestep()
error += qresidual[i] * qresidual[i];
// on the first iteration, compute initial error
- if (numIter == 0)
+ if (numIter == 0)
{
error0 = error;
errorQuotient = 1.0;
@@ -289,7 +289,7 @@ int ImplicitNewmarkSparse::DoTimestep()
else
{
// error divided by the initial error, before performing this iteration
- errorQuotient = error / error0;
+ errorQuotient = error / error0;
}
if (errorQuotient < epsilon * epsilon)
@@ -381,16 +381,16 @@ int ImplicitNewmarkSparse::DoTimestep()
}
void ImplicitNewmarkSparse::UseStaticSolver(bool useStaticSolver_)
-{
+{
useStaticSolver = useStaticSolver_;
- if (!useStaticSolver)
+ if (!useStaticSolver)
{
memset(qvel, 0, sizeof(double) * r);
memset(qaccel, 0, sizeof(double) * r);
- memset(qvel_1, 0, sizeof(double) * r);
+ memset(qvel_1, 0, sizeof(double) * r);
memset(qaccel_1, 0, sizeof(double) * r);
memcpy(q_1, q, sizeof(double) * r);
}
-}
+}
diff --git a/src/libintegratorSparse/implicitNewmarkSparse.h b/src/libintegratorSparse/implicitNewmarkSparse.h
index 301793b3cee06ba8112fc4282eb80c8b17f4db0f..095d574ae44985627d53f62c88f5f513c9cb8721 100644
--- a/src/libintegratorSparse/implicitNewmarkSparse.h
+++ b/src/libintegratorSparse/implicitNewmarkSparse.h
@@ -32,7 +32,7 @@
See also integratorBase.h .
- This class either uses SPOOLES, PARDISO, or our own Jacobi-preconitioned
+ This class either uses SPOOLES, PARDISO, or our own Jacobi-preconitioned
CG to solve the large sparse linear systems.
You can switch between these solvers at compile time,
@@ -76,7 +76,7 @@ public:
// constrainedDOFs is an integer array of degrees of freedom that are to be fixed to zero (e.g., to permanently fix a vertex in a deformable simulation)
// constrainedDOFs are 0-indexed (separate DOFs for x,y,z), and must be pre-sorted (ascending)
// numThreads applies only to the PARDISO solver; if numThreads > 0, the sparse linear solves are multi-threaded; default: 0 (use single-threading)
- ImplicitNewmarkSparse(int r, double timestep, SparseMatrix * massMatrix, ForceModel * forceModel, int positiveDefiniteSolver=0, int numConstrainedDOFs=0, int * constrainedDOFs=NULL, double dampingMassCoef=0.0, double dampingStiffnessCoef=0.0, int maxIterations = 1, double epsilon = 1E-6, double NewmarkBeta=0.25, double NewmarkGamma=0.5, int numSolverThreads=0);
+ ImplicitNewmarkSparse(int r, double timestep, SparseMatrix * massMatrix, ForceModel * forceModel, int positiveDefiniteSolver=0, int numConstrainedDOFs=0, int * constrainedDOFs=NULL, double dampingMassCoef=0.0, double dampingStiffnessCoef=0.0, int maxIterations = 1, double epsilon = 1E-6, double NewmarkBeta=0.25, double NewmarkGamma=0.5, int numSolverThreads=0);
virtual ~ImplicitNewmarkSparse();
@@ -84,7 +84,7 @@ public:
virtual void SetDampingMatrix(SparseMatrix * dampingMatrix);
inline virtual void SetTimestep(double timestep) { this->timestep = timestep; UpdateAlphas(); }
- // sets q, qvel
+ // sets q, qvel
// automatically computes acceleration assuming zero external force
// returns 0 on succes, 1 if solver fails to converge
// note: there are also other state setting routines in the base class
@@ -92,7 +92,7 @@ public:
// performs one step of simulation (returns 0 on sucess, and 1 on failure)
// failure can occur, for example, if you are using the positive definite solver and the system matrix has negative eigenvalues
- virtual int DoTimestep();
+ virtual int DoTimestep();
inline void SetNewmarkBeta(double NewmarkBeta) { this->NewmarkBeta = NewmarkBeta; UpdateAlphas(); }
inline void SetNewmarkGamma(double NewmarkGamma) { this->NewmarkGamma = NewmarkGamma; UpdateAlphas(); }
@@ -110,7 +110,7 @@ protected:
// parameters for implicit Newmark
double NewmarkBeta,NewmarkGamma;
double alpha1, alpha2, alpha3, alpha4, alpha5, alpha6;
- double epsilon;
+ double epsilon;
int maxIterations;
void UpdateAlphas();
diff --git a/src/libintegratorSparse/integratorBaseSparse.h b/src/libintegratorSparse/integratorBaseSparse.h
index 56eedb907c7444ae33b2632a7c76ba5ac0b372c1..2d72119f49fe3069a94996bb41acd0b957c7a20e 100644
--- a/src/libintegratorSparse/integratorBaseSparse.h
+++ b/src/libintegratorSparse/integratorBaseSparse.h
@@ -68,7 +68,7 @@ public:
virtual double GetTotalMass();
protected:
- SparseMatrix * massMatrix;
+ SparseMatrix * massMatrix;
ForceModel * forceModel;
int ownDampingMatrix;
SparseMatrix * dampingMatrix;
diff --git a/src/libsparseSolver/ARPACKSolver.cpp b/src/libsparseSolver/ARPACKSolver.cpp
index c07491cd56d23991ae465b704666a41b4c9ab918..fdd1a5b28f61cfa35fecd2f3d7954336ed783f4f 100644
--- a/src/libsparseSolver/ARPACKSolver.cpp
+++ b/src/libsparseSolver/ARPACKSolver.cpp
@@ -48,7 +48,7 @@ using namespace std;
#define DSEUPD dseupd_
#define INTEGER int
-/* Note:
+/* Note:
This code is now independent of ARPACK++ and calls ARPACK directly. It has been tested and works fine.
Previously, it used ARPACK++, which, while convenient, required installing ARPACK++ (which required
@@ -68,7 +68,7 @@ using namespace std;
double * TOL, double * RESID, INTEGER * NCV, double * V, INTEGER * LDV, INTEGER * IPARAM,
INTEGER * IPNTR, double * WORKD, double * WORKL, INTEGER * LWORKL, INTEGER * INFO);
- // call DSAUPD
+ // call DSAUPD
//c ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM,
//c IPNTR, WORKD, WORKL, LWORKL, INFO )
@@ -78,7 +78,7 @@ using namespace std;
double * TOL, double * RESID, INTEGER * NCV, double * V, INTEGER * LDV, INTEGER * IPARAM,
INTEGER * IPNTR, double * WORKD, double * WORKL, INTEGER * LWORKL, INTEGER * INFO);
- //c call DSEUPD
+ //c call DSEUPD
//c ( RVEC, HOWMNY, SELECT, D, Z, LDZ, SIGMA, BMAT, N, WHICH, NEV, TOL,
//c RESID, NCV, V, LDV, IPARAM, IPNTR, WORKD, WORKL, LWORKL, INFO )
}
@@ -140,7 +140,7 @@ int ARPACKSolver::SolveGenEigReg(SparseMatrix * KInput, SparseMatrix * MInput, i
#ifdef PARDISO
if (verbose >= 1)
printf("Linear solver: PARDISO (%d threads).\n", (numLinearSolverThreads == 0) ? 1 : numLinearSolverThreads);
- int positiveDefinite = 0;
+ int positiveDefinite = 0;
int directIterative = 0;
PardisoSolver * pardisoSolver = new PardisoSolver(M, numLinearSolverThreads, positiveDefinite, directIterative, verbose);
pardisoSolver->ComputeCholeskyDecomposition(M);
@@ -165,7 +165,7 @@ int ARPACKSolver::SolveGenEigReg(SparseMatrix * KInput, SparseMatrix * MInput, i
int maxitp = 0, ARFLOAT* residp = NULL, bool ishiftp = true);
*/
ARSymGenEig<double, InvMKSolver, SparseMatrix> solver
- (np, numEigenvalues,
+ (np, numEigenvalues,
&invMKSolver, &InvMKSolver::ComputeInvMK,
M, (void (SparseMatrix::*)(double*,double*)) &SparseMatrix::MultiplyVector, "LM");
@@ -177,7 +177,7 @@ int ARPACKSolver::SolveGenEigReg(SparseMatrix * KInput, SparseMatrix * MInput, i
INTEGER IDO = 0;
char BMAT = 'G';
INTEGER N = np;
- char WHICH[3] = "LM";
+ char WHICH[3] = "LM";
INTEGER NEV = numEigenvalues;
double TOL = 0.0;
double * RESID = (double*) malloc (sizeof(double) * N);
@@ -188,18 +188,18 @@ int ARPACKSolver::SolveGenEigReg(SparseMatrix * KInput, SparseMatrix * MInput, i
INTEGER LDV = N;
INTEGER IPARAM[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
IPARAM[0] = 1; // no user-provided shifts
- IPARAM[2] = maxIter;
+ IPARAM[2] = maxIter;
IPARAM[3] = 1;
//IPARAM[4] = 0;
IPARAM[6] = 2; // the mode (regular generalized eigenproblem)
- INTEGER IPNTR[11];
+ INTEGER IPNTR[11];
double * WORKD = (double*) malloc (sizeof(double) * 3 * N);
- INTEGER LWORKL = 2 * (NCV*NCV + 8*NCV); // at least NCV**2 + 8*NCV;
+ INTEGER LWORKL = 2 * (NCV*NCV + 8*NCV); // at least NCV**2 + 8*NCV;
double * WORKL = (double*) malloc (sizeof(double) * LWORKL);
INTEGER INFO = 0;
double * buffer = (double*) malloc (sizeof(double) * N);
-
+
do
{
//printf("Entering dsaupd...\n"); fflush(NULL);
@@ -212,7 +212,7 @@ int ARPACKSolver::SolveGenEigReg(SparseMatrix * KInput, SparseMatrix * MInput, i
}
//printf("IDO=%d\n", (int)IDO); fflush(NULL);
-
+
switch(IDO)
{
case -1:
@@ -255,8 +255,8 @@ int ARPACKSolver::SolveGenEigReg(SparseMatrix * KInput, SparseMatrix * MInput, i
double * SIGMA = NULL;
DSEUPD(&RVEC, &HOWMNY, SELECT, D, Z, &LDZ, SIGMA,
&BMAT, &N, WHICH, &NEV, &TOL, RESID, &NCV, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO);
-
- free(SELECT);
+
+ free(SELECT);
free(RESID);
free(V);
free(WORKD);
@@ -315,8 +315,8 @@ int ARPACKSolver::SolveGenEigReg(SparseMatrix * KInput, SparseMatrix * MInput, i
// Printing eigenvalues.
cout << "Eigenvalues:" << endl;
- //for (int i=numEigenvalues-nconv; i<numEigenvalues; i++)
- for (int i=0; i<nconv; i++)
+ //for (int i=numEigenvalues-nconv; i<numEigenvalues; i++)
+ for (int i=0; i<nconv; i++)
{
cout << " lambda[" << (i+1) << "]: " << eigenvalues[i] << endl;
}
@@ -335,7 +335,7 @@ int ARPACKSolver::SolveGenEigReg(SparseMatrix * KInput, SparseMatrix * MInput, i
if (infinityNormM > infNorm)
infNorm = infinityNormM;
- for (int i=0; i<nconv; i++)
+ for (int i=0; i<nconv; i++)
{
KInput->MultiplyVector(&eigenvectors[np*i],Ax);
MInput->MultiplyVector(&eigenvectors[np*i],Bx);
@@ -350,7 +350,7 @@ int ARPACKSolver::SolveGenEigReg(SparseMatrix * KInput, SparseMatrix * MInput, i
}
/*
- for (int i=0; i<nconv; i++)
+ for (int i=0; i<nconv; i++)
{
cout << "||A*x(" << (i+1) << ") - lambda(" << (i+1);
cout << ")*B*x(" << (i+1) << ")|| / |lambda| / max(||A||,||B||) = " << ResNorm[i] << endl;
@@ -382,7 +382,7 @@ int ARPACKSolver::SolveGenEigShInv(SparseMatrix * K, SparseMatrix * M, int numEi
if (sigma != 0)
{
// compute KsigmaM = K - sigma * M
- KsigmaM = new SparseMatrix(*K);
+ KsigmaM = new SparseMatrix(*K);
KsigmaM->BuildSubMatrixIndices(*M);
KsigmaM->AddSubMatrix(-sigma, *M);
}
@@ -408,7 +408,7 @@ int ARPACKSolver::SolveGenEigShInv(SparseMatrix * K, SparseMatrix * M, int numEi
#ifdef PARDISO
if (verbose >= 1)
printf("Linear solver: PARDISO (%d threads).\n", (numLinearSolverThreads == 0) ? 1 : numLinearSolverThreads);
- int positiveDefinite = 0;
+ int positiveDefinite = 0;
int directIterative = 0;
PardisoSolver * pardisoSolver = new PardisoSolver(KsigmaM, numLinearSolverThreads, positiveDefinite, directIterative, verbose);
pardisoSolver->ComputeCholeskyDecomposition(KsigmaM);
@@ -430,7 +430,7 @@ int ARPACKSolver::SolveGenEigShInv(SparseMatrix * K, SparseMatrix * M, int numEi
INTEGER IDO = 0;
char BMAT = 'G';
INTEGER N = np;
- char WHICH[3] = "LM";
+ char WHICH[3] = "LM";
INTEGER NEV = numEigenvalues;
double TOL = 0.0;
double * RESID = (double*) malloc (sizeof(double) * N);
@@ -441,17 +441,17 @@ int ARPACKSolver::SolveGenEigShInv(SparseMatrix * K, SparseMatrix * M, int numEi
INTEGER LDV = N;
INTEGER IPARAM[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
IPARAM[0] = 1; // no user-provided shifts
- IPARAM[2] = maxIter;
+ IPARAM[2] = maxIter;
IPARAM[3] = 1;
IPARAM[6] = 3; // the mode (shift-inverted generalized eigenproblem)
- INTEGER IPNTR[11];
+ INTEGER IPNTR[11];
double * WORKD = (double*) malloc (sizeof(double) * 3 * N);
- INTEGER LWORKL = 2 * (NCV*NCV + 8*NCV); // at least NCV**2 + 8*NCV;
+ INTEGER LWORKL = 2 * (NCV*NCV + 8*NCV); // at least NCV**2 + 8*NCV;
double * WORKL = (double*) malloc (sizeof(double) * LWORKL);
INTEGER INFO = 0;
double * buffer = (double*) malloc (sizeof(double) * N);
-
+
if (verbose >= 1)
{
cout << "Dimension of the system : " << np << endl;
@@ -471,8 +471,8 @@ int ARPACKSolver::SolveGenEigShInv(SparseMatrix * K, SparseMatrix * M, int numEi
}
//printf("IDO=%d\n", (int)IDO); fflush(NULL);
- //c ===> OP = (inv[K - sigma*M])*M and B = M.
-
+ //c ===> OP = (inv[K - sigma*M])*M and B = M.
+
switch(IDO)
{
case -1:
@@ -530,8 +530,8 @@ int ARPACKSolver::SolveGenEigShInv(SparseMatrix * K, SparseMatrix * M, int numEi
double SIGMA = sigma;
DSEUPD(&RVEC, &HOWMNY, SELECT, D, Z, &LDZ, &SIGMA,
&BMAT, &N, WHICH, &NEV, &TOL, RESID, &NCV, V, &LDV, IPARAM, IPNTR, WORKD, WORKL, &LWORKL, &INFO);
-
- free(SELECT);
+
+ free(SELECT);
free(RESID);
free(V);
free(WORKD);
@@ -565,8 +565,8 @@ int ARPACKSolver::SolveGenEigShInv(SparseMatrix * K, SparseMatrix * M, int numEi
// Printing eigenvalues.
cout << "Eigenvalues:" << endl;
- //for (int i=numEigenvalues-nconv; i<numEigenvalues; i++)
- for (int i=0; i<nconv; i++)
+ //for (int i=numEigenvalues-nconv; i<numEigenvalues; i++)
+ for (int i=0; i<nconv; i++)
{
cout << " lambda[" << (i+1) << "]: " << eigenvalues[i] << endl;
}
@@ -585,7 +585,7 @@ int ARPACKSolver::SolveGenEigShInv(SparseMatrix * K, SparseMatrix * M, int numEi
if (infinityNormM > infNorm)
infNorm = infinityNormM;
- for (int i=0; i<nconv; i++)
+ for (int i=0; i<nconv; i++)
{
K->MultiplyVector(&eigenvectors[np*i],Ax);
M->MultiplyVector(&eigenvectors[np*i],Bx);
@@ -600,7 +600,7 @@ int ARPACKSolver::SolveGenEigShInv(SparseMatrix * K, SparseMatrix * M, int numEi
}
/*
- for (int i=0; i<nconv; i++)
+ for (int i=0; i<nconv; i++)
{
cout << "||A*x(" << (i+1) << ") - lambda(" << (i+1);
cout << ")*B*x(" << (i+1) << ")|| / |lambda| / max(||A||,||B||) = " << ResNorm[i] << endl;
@@ -626,31 +626,31 @@ c\BeginDoc
c
c\Name: dsaupd
c
-c\Description:
+c\Description:
c
-c Reverse communication interface for the Implicitly Restarted Arnoldi
-c Iteration. For symmetric problems this reduces to a variant of the Lanczos
-c method. This method has been designed to compute approximations to a
-c few eigenpairs of a linear operator OP that is real and symmetric
-c with respect to a real positive semi-definite symmetric matrix B,
+c Reverse communication interface for the Implicitly Restarted Arnoldi
+c Iteration. For symmetric problems this reduces to a variant of the Lanczos
+c method. This method has been designed to compute approximations to a
+c few eigenpairs of a linear operator OP that is real and symmetric
+c with respect to a real positive semi-definite symmetric matrix B,
c i.e.
-c
-c B*OP = (OP')*B.
c
-c Another way to express this condition is
+c B*OP = (OP')*B.
+c
+c Another way to express this condition is
c
c < x,OPy > = < OPx,y > where < z,w > = z'Bw .
-c
-c In the standard eigenproblem B is the identity matrix.
+c
+c In the standard eigenproblem B is the identity matrix.
c ( A' denotes transpose of A)
c
c The computed approximate eigenvalues are called Ritz values and
c the corresponding approximate eigenvectors are called Ritz vectors.
c
-c dsaupd is usually called iteratively to solve one of the
+c dsaupd is usually called iteratively to solve one of the
c following problems:
c
-c Mode 1: A*x = lambda*x, A symmetric
+c Mode 1: A*x = lambda*x, A symmetric
c ===> OP = A and B = I.
c
c Mode 2: A*x = lambda*M*x, A symmetric, M symmetric positive definite
@@ -658,10 +658,10 @@ c ===> OP = inv[M]*A and B = M.
c ===> (If M can be factored see remark 3 below)
c
c Mode 3: K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite
-c ===> OP = (inv[K - sigma*M])*M and B = M.
+c ===> OP = (inv[K - sigma*M])*M and B = M.
c ===> Shift-and-Invert mode
c
-c Mode 4: K*x = lambda*KG*x, K symmetric positive semi-definite,
+c Mode 4: K*x = lambda*KG*x, K symmetric positive semi-definite,
c KG symmetric indefinite
c ===> OP = (inv[K - sigma*KG])*K and B = K.
c ===> Buckling mode
@@ -683,13 +683,13 @@ c the accuracy requirements for the eigenvalue
c approximations.
c
c\Usage:
-c call dsaupd
+c call dsaupd
c ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM,
c IPNTR, WORKD, WORKL, LWORKL, INFO )
c
c\Arguments
c IDO Integer. (INPUT/OUTPUT)
-c Reverse communication flag. IDO must be zero on the first
+c Reverse communication flag. IDO must be zero on the first
c call to dsaupd. IDO will be set internally to
c indicate the type of operation to be performed. Control is
c then given back to the calling routine which has the
@@ -716,11 +716,11 @@ c IPNTR(11) is the pointer into WORKL for
c placing the shifts. See remark 6 below.
c IDO = 99: done
c -------------------------------------------------------------
-c After the initialization phase, when the routine is used in
+c After the initialization phase, when the routine is used in
c either the "shift-and-invert" mode or the Cayley transform
-c mode, the vector B * X is already available and does not
+c mode, the vector B * X is already available and does not
c need to be recomputed in forming OP*X.
-c
+c
c BMAT Character*1. (INPUT)
c BMAT specifies the type of the matrix B that defines the
c semi-inner product for the operator OP.
@@ -736,7 +736,7 @@ c
c 'LA' - compute the NEV largest (algebraic) eigenvalues.
c 'SA' - compute the NEV smallest (algebraic) eigenvalues.
c 'LM' - compute the NEV largest (in magnitude) eigenvalues.
-c 'SM' - compute the NEV smallest (in magnitude) eigenvalues.
+c 'SM' - compute the NEV smallest (in magnitude) eigenvalues.
c 'BE' - compute NEV eigenvalues, half from each end of the
c spectrum. When NEV is odd, compute one more from the
c high end than from the low end.
@@ -746,27 +746,27 @@ c NEV Integer. (INPUT)
c Number of eigenvalues of OP to be computed. 0 < NEV < N.
c
c TOL Double precision scalar. (INPUT)
-c Stopping criterion: the relative accuracy of the Ritz value
+c Stopping criterion: the relative accuracy of the Ritz value
c is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I)).
c If TOL .LE. 0. is passed a default is set:
c DEFAULT = DLAMCH('EPS') (machine precision as computed
c by the LAPACK auxiliary subroutine DLAMCH).
c
c RESID Double precision array of length N. (INPUT/OUTPUT)
-c On INPUT:
+c On INPUT:
c If INFO .EQ. 0, a random initial residual vector is used.
c If INFO .NE. 0, RESID contains the initial residual vector,
c possibly from a previous run.
c On OUTPUT:
-c RESID contains the final residual vector.
+c RESID contains the final residual vector.
c
c NCV Integer. (INPUT)
c Number of columns of the matrix V (less than or equal to N).
-c This will indicate how many Lanczos vectors are generated
-c at each iteration. After the startup phase in which NEV
-c Lanczos vectors are generated, the algorithm generates
+c This will indicate how many Lanczos vectors are generated
+c at each iteration. After the startup phase in which NEV
+c Lanczos vectors are generated, the algorithm generates
c NCV-NEV Lanczos vectors at each subsequent update iteration.
-c Most of the cost in generating each Lanczos vector is in the
+c Most of the cost in generating each Lanczos vector is in the
c matrix-vector product OP*x. (See remark 4 below).
c
c V Double precision N by NCV array. (OUTPUT)
@@ -786,10 +786,10 @@ c reverse communication. The NCV eigenvalues of
c the current tridiagonal matrix T are returned in
c the part of WORKL array corresponding to RITZ.
c See remark 6 below.
-c ISHIFT = 1: exact shifts with respect to the reduced
-c tridiagonal matrix T. This is equivalent to
-c restarting the iteration with a starting vector
-c that is a linear combination of Ritz vectors
+c ISHIFT = 1: exact shifts with respect to the reduced
+c tridiagonal matrix T. This is equivalent to
+c restarting the iteration with a starting vector
+c that is a linear combination of Ritz vectors
c associated with the "wanted" Ritz values.
c -------------------------------------------------------------
c
@@ -797,8 +797,8 @@ c IPARAM(2) = LEVEC
c No longer referenced. See remark 2 below.
c
c IPARAM(3) = MXITER
-c On INPUT: maximum number of Arnoldi update iterations allowed.
-c On OUTPUT: actual number of Arnoldi update iterations taken.
+c On INPUT: maximum number of Arnoldi update iterations allowed.
+c On OUTPUT: actual number of Arnoldi update iterations taken.
c
c IPARAM(4) = NB: blocksize to be used in the recurrence.
c The code currently works only for NB = 1.
@@ -808,11 +808,11 @@ c This represents the number of Ritz values that satisfy
c the convergence criterion.
c
c IPARAM(6) = IUPD
-c No longer referenced. Implicit restarting is ALWAYS used.
+c No longer referenced. Implicit restarting is ALWAYS used.
c
c IPARAM(7) = MODE
c On INPUT determines what type of eigenproblem is being solved.
-c Must be 1,2,3,4,5; See under \Description of dsaupd for the
+c Must be 1,2,3,4,5; See under \Description of dsaupd for the
c five modes available.
c
c IPARAM(8) = NP
@@ -824,7 +824,7 @@ c
c IPARAM(9) = NUMOP, IPARAM(10) = NUMOPB, IPARAM(11) = NUMREO,
c OUTPUT: NUMOP = total number of OP*x operations,
c NUMOPB = total number of B*x operations if BMAT='G',
-c NUMREO = total number of steps of re-orthogonalization.
+c NUMREO = total number of steps of re-orthogonalization.
c
c IPNTR Integer array of length 11. (OUTPUT)
c Pointer to mark the starting locations in the WORKD and WORKL
@@ -832,7 +832,7 @@ c arrays for matrices/vectors used by the Lanczos iteration.
c -------------------------------------------------------------
c IPNTR(1): pointer to the current operand vector X in WORKD.
c IPNTR(2): pointer to the current result vector Y in WORKD.
-c IPNTR(3): pointer to the vector B * X in WORKD when used in
+c IPNTR(3): pointer to the vector B * X in WORKD when used in
c the shift-and-invert mode.
c IPNTR(4): pointer to the next available location in WORKL
c that is untouched by the program.
@@ -849,14 +849,14 @@ c dseupd if RVEC = .TRUE. See Remarks.
c Note: IPNTR(8:10) is only referenced by dseupd. See Remark 2.
c IPNTR(11): pointer to the NP shifts in WORKL. See Remark 6 below.
c -------------------------------------------------------------
-c
+c
c WORKD Double precision work array of length 3*N. (REVERSE COMMUNICATION)
c Distributed array to be used in the basic Arnoldi iteration
-c for reverse communication. The user should not use WORKD
+c for reverse communication. The user should not use WORKD
c as temporary workspace during the iteration. Upon termination
c WORKD(1:N) contains B*RESID(1:N). If the Ritz vectors are desired
c subroutine dseupd uses this output.
-c See Data Distribution Note below.
+c See Data Distribution Note below.
c
c WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE)
c Private (replicated) array on each PE or array allocated on
@@ -872,13 +872,13 @@ c possibly from a previous run.
c Error flag on output.
c = 0: Normal exit.
c = 1: Maximum number of iterations taken.
-c All possible eigenvalues of OP has been found. IPARAM(5)
+c All possible eigenvalues of OP has been found. IPARAM(5)
c returns the number of wanted converged Ritz values.
c = 2: No longer an informational error. Deprecated starting
c with release 2 of ARPACK.
-c = 3: No shifts could be applied during a cycle of the
-c Implicitly restarted Arnoldi iteration. One possibility
-c is to increase the size of NCV relative to NEV.
+c = 3: No shifts could be applied during a cycle of the
+c Implicitly restarted Arnoldi iteration. One possibility
+c is to increase the size of NCV relative to NEV.
c See remark 4 below.
c = -1: N must be positive.
c = -2: NEV must be positive.
@@ -902,12 +902,12 @@ c enough workspace and array storage has been allocated.
c
c
c\Remarks
-c 1. The converged Ritz values are always returned in ascending
+c 1. The converged Ritz values are always returned in ascending
c algebraic order. The computed Ritz values are approximate
c eigenvalues of OP. The selection of WHICH should be made
-c with this in mind when Mode = 3,4,5. After convergence,
-c approximate eigenvalues of the original problem may be obtained
-c with the ARPACK subroutine dseupd.
+c with this in mind when Mode = 3,4,5. After convergence,
+c approximate eigenvalues of the original problem may be obtained
+c with the ARPACK subroutine dseupd.
c
c 2. If the Ritz vectors corresponding to the converged Ritz values
c are needed, the user must call dseupd immediately following completion
@@ -915,7 +915,7 @@ c of dsaupd. This is new starting with version 2.1 of ARPACK.
c
c 3. If M can be factored into a Cholesky factorization M = LL'
c then Mode = 2 should not be selected. Instead one should use
-c Mode = 1 with OP = inv(L)*A*inv(L'). Appropriate triangular
+c Mode = 1 with OP = inv(L)*A*inv(L'). Appropriate triangular
c linear systems should be solved with L and L' rather
c than computing inverses. After convergence, an approximate
c eigenvector z of the original problem is recovered by solving
@@ -925,7 +925,7 @@ c 4. At present there is no a-priori analysis to guide the selection
c of NCV relative to NEV. The only formal requirement is that NCV > NEV.
c However, it is recommended that NCV .ge. 2*NEV. If many problems of
c the same type are to be solved, one should experiment with increasing
-c NCV while keeping NEV fixed for a given test problem. This will
+c NCV while keeping NEV fixed for a given test problem. This will
c usually decrease the required number of OP*x operations but it
c also increases the work and storage required to maintain the orthogonal
c basis vectors. The optimal "cross-over" with respect to CPU time
@@ -937,16 +937,16 @@ c When IPARAM(7) = 2 OP = inv(B)*A. After computing A*X the user
c must overwrite X with A*X. Y is then the solution to the linear set
c of equations B*Y = A*X.
c
-c 6. When IPARAM(1) = 0, and IDO = 3, the user needs to provide the
-c NP = IPARAM(8) shifts in locations:
-c 1 WORKL(IPNTR(11))
-c 2 WORKL(IPNTR(11)+1)
-c .
-c .
-c .
-c NP WORKL(IPNTR(11)+NP-1).
+c 6. When IPARAM(1) = 0, and IDO = 3, the user needs to provide the
+c NP = IPARAM(8) shifts in locations:
+c 1 WORKL(IPNTR(11))
+c 2 WORKL(IPNTR(11)+1)
+c .
+c .
+c .
+c NP WORKL(IPNTR(11)+NP-1).
c
-c The eigenvalues of the current tridiagonal matrix are located in
+c The eigenvalues of the current tridiagonal matrix are located in
c WORKL(IPNTR(6)) through WORKL(IPNTR(6)+NCV-1). They are in the
c order defined by WHICH. The associated Ritz estimates are located in
c WORKL(IPNTR(8)), WORKL(IPNTR(8)+1), ... , WORKL(IPNTR(8)+NCV-1).
diff --git a/src/libsparseSolver/ARPACKSolver.h b/src/libsparseSolver/ARPACKSolver.h
index 8531095f1cc5b92f6da3b52d9069ba300fa994c8..ed61d375fd7fb0552b6078ae15c1b59e8e0d3d42 100644
--- a/src/libsparseSolver/ARPACKSolver.h
+++ b/src/libsparseSolver/ARPACKSolver.h
@@ -39,7 +39,7 @@
class ARPACKSolver
{
public:
-
+
// K * x = lambda * M * x
// returns the number of converged eigenvalues
// assumes that both K and M are symmetric, and that M > 0
diff --git a/src/libsparseSolver/CGSolver.cpp b/src/libsparseSolver/CGSolver.cpp
index f43b93d681375c651de9d75f941f4c3935cc51ac..8a1fb846ca965262c56a7ad8d2205848b11914c8 100644
--- a/src/libsparseSolver/CGSolver.cpp
+++ b/src/libsparseSolver/CGSolver.cpp
@@ -31,7 +31,7 @@
#include <stdio.h>
#include "CGSolver.h"
-CGSolver::CGSolver(SparseMatrix * A_): A(A_)
+CGSolver::CGSolver(SparseMatrix * A_): A(A_)
{
numRows = A->GetNumRows();
InitBuffers();
@@ -78,7 +78,7 @@ void CGSolver::InitBuffers()
}
// implements the virtual method from LinearSolver by calling "SolveLinearSystem" with default parameters
-int CGSolver::SolveLinearSystem(double * x, const double * b)
+int CGSolver::SolveLinearSystem(double * x, const double * b)
{
return SolveLinearSystemWithJacobiPreconditioner(x, b, 1E-6, 1000, 0);
}
diff --git a/src/libsparseSolver/CGSolver.h b/src/libsparseSolver/CGSolver.h
index 8ff9718cda312d4bcf83c5a48fb3b35c0c797fef..70f4ee6591803109d28c105c5a48f7d1583ac6e2 100755
--- a/src/libsparseSolver/CGSolver.h
+++ b/src/libsparseSolver/CGSolver.h
@@ -28,7 +28,7 @@
/*
A conjugate gradient solver built on top of the sparse matrix class.
- There are two solver versions: without preconditioning, and with
+ There are two solver versions: without preconditioning, and with
Jacobi preconditioning.
You can either provide a sparse matrix, or a callback function to
@@ -36,7 +36,7 @@
The sparse matrix must be symmetric and positive-definite.
- The CG solvers were implemented by following Jonathan Shewchuk's
+ The CG solvers were implemented by following Jonathan Shewchuk's
An Introduction to the Conjugate Gradient Method Without the Agonizing Pain:
http://www.cs.cmu.edu/~jrs/jrspapers.html#cg
@@ -57,7 +57,7 @@ public:
CGSolver(SparseMatrix * A);
// This constructor makes it possible to only provide a
- // "black-box" matrix-vector multiplication routine
+ // "black-box" matrix-vector multiplication routine
// (no need to explicitly give the matrix):
// given x, the routine must compute A * x, and store it into Ax.
// "data" should not be used/touched by the user-written "black-box" routine.
@@ -74,7 +74,7 @@ public:
// output: solution (in x)
// "eps" is the convergence criterium: solver converges when the L2 residual errors is less than eps times the initial L2 residual error, must have 0 < eps < 1
// maximum number of conjugate-gradient iterations is set by "maxIterations"
- // return value is the number of iterations performed
+ // return value is the number of iterations performed
// if solver did not converge, the return value will have a negative sign
int SolveLinearSystemWithoutPreconditioner(double * x, const double * b, double eps=1e-6, int maxIterations=1000, int verbose=0);
@@ -91,7 +91,7 @@ protected:
int numRows;
blackBoxProductType multiplicator;
void * multiplicatorData;
- SparseMatrix * A;
+ SparseMatrix * A;
double * r, * d, * q; // terminology from Shewchuk's work
double * invDiagonal;
diff --git a/src/libsparseSolver/PardisoSolver.cpp b/src/libsparseSolver/PardisoSolver.cpp
index ac5c4c5a920e3c074a350124d342ebcd08bac00f..02cae2b62c03cfcbdc292b2cebfc9678f567acf0 100644
--- a/src/libsparseSolver/PardisoSolver.cpp
+++ b/src/libsparseSolver/PardisoSolver.cpp
@@ -70,9 +70,9 @@ PardisoSolver::PardisoSolver(const SparseMatrix * A, int numThreads_, int positi
int numUpperTriangleEntries = A->GetNumUpperTriangleEntries();
- a = (double*) malloc (sizeof(double) * numUpperTriangleEntries);
- ia = (int*) malloc (sizeof(int) * (A->GetNumRows() + 1));
- ja = (int*) malloc (sizeof(int) * numUpperTriangleEntries);
+ a = (double*) malloc (sizeof(double) * numUpperTriangleEntries);
+ ia = (int*) malloc (sizeof(int) * (A->GetNumRows() + 1));
+ ja = (int*) malloc (sizeof(int) * numUpperTriangleEntries);
int upperTriangleOnly = 1;
int oneIndexed = 1;
@@ -90,7 +90,7 @@ PardisoSolver::PardisoSolver(const SparseMatrix * A, int numThreads_, int positi
msglvl = verbose >= 1 ? verbose-1 : 0; /* Print statistical information in file */
error = 0; /* Initialize error flag */
- for (int i = 0; i < 64; i++)
+ for (int i = 0; i < 64; i++)
iparm[i] = 0;
iparm[0] = 1; // No solver default
iparm[1] = 2; // 0=minimum degree ordering, 2=Fill-in reordering from METIS
@@ -123,7 +123,7 @@ PardisoSolver::PardisoSolver(const SparseMatrix * A, int numThreads_, int positi
/* .. Initialize the internal solver memory pointer. This is only */
/* necessary for the FIRST call of the PARDISO solver. */
/* -------------------------------------------------------------------- */
- for (int i=0; i<64; i++)
+ for (int i=0; i<64; i++)
pt[i] = 0;
if (verbose >= 1)
@@ -180,7 +180,7 @@ MKL_INT PardisoSolver::ComputeCholeskyDecomposition(const SparseMatrix * A)
int oneIndexed = 1;
A->GenerateCompressedRowMajorFormat(a, NULL, NULL, upperTriangleOnly, oneIndexed);
- // factor
+ // factor
phase = 22;
PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, NULL, &nrhs, iparm, &msglvl, NULL, NULL, &error);
@@ -211,7 +211,7 @@ int PardisoSolver::SolveLinearSystem(double * x, const double * rhs)
printf("Error: Pardiso solve returned non-zero exit code %d.\n", error);
if (verbose >= 2)
- printf("Solve completed.\n");
+ printf("Solve completed.\n");
return (int)error;
}
@@ -240,7 +240,7 @@ MKL_INT PardisoSolver::ForwardSubstitution(double * y, const double * rhs)
printf("Error: Pardiso solve returned non-zero exit code %d.\n", error);
if (verbose >= 2)
- printf("Solve completed.\n");
+ printf("Solve completed.\n");
return error;
}
@@ -274,7 +274,7 @@ MKL_INT PardisoSolver::DiagonalSubstitution(double * v, const double * y)
printf("Error: Pardiso solve returned non-zero exit code %d.\n", error);
if (verbose >= 2)
- printf("Solve completed.\n");
+ printf("Solve completed.\n");
return error;
}
@@ -302,7 +302,7 @@ MKL_INT PardisoSolver::BackwardSubstitution(double * x, const double * y)
printf("Error: Pardiso solve returned non-zero exit code %d.\n", error);
if (verbose >= 2)
- printf("Solve completed.\n");
+ printf("Solve completed.\n");
return error;
}
@@ -353,7 +353,7 @@ MKL_INT PardisoSolver::SolveLinearSystemDirectIterative(const SparseMatrix * A,
printf("Error: Pardiso solve returned non-zero exit code %d.\n", error);
if (verbose >= 2)
- printf("Solve completed.\n");
+ printf("Solve completed.\n");
return error;
}
@@ -362,7 +362,7 @@ MKL_INT PardisoSolver::SolveLinearSystemDirectIterative(const SparseMatrix * A,
// Pardiso Solver is not available
-PardisoSolver::PardisoSolver(const SparseMatrix * A, int numThreads_, int positiveDefinite_, int directIterative_, int verbose_) : numThreads(numThreads_), positiveDefinite(positiveDefinite_), directIterative(directIterative_), verbose(verbose_)
+PardisoSolver::PardisoSolver(const SparseMatrix * A, int numThreads_, int positiveDefinite_, int directIterative_, int verbose_) : numThreads(numThreads_), positiveDefinite(positiveDefinite_), directIterative(directIterative_), verbose(verbose_)
{
DisabledSolverError();
}
diff --git a/src/libsparseSolver/PardisoSolver.h b/src/libsparseSolver/PardisoSolver.h
index be13bda0ffaafaa5261c5a21db06686d403d2488..e9d97490789f04b3bea3c9c0a03ab1e8768bcfc8 100644
--- a/src/libsparseSolver/PardisoSolver.h
+++ b/src/libsparseSolver/PardisoSolver.h
@@ -35,7 +35,7 @@
Uses a multi-threaded approach to perform the solve.
The solution is obtained using the the Pardiso library .
-
+
Jernej Barbic, MIT, 2007-2009
*/
@@ -49,7 +49,7 @@ class PardisoSolver : public LinearSolver
{
public:
- // the constructor will compute the permutation to re-order A
+ // the constructor will compute the permutation to re-order A
// only the topology of A matters for this step
// A is not modified
PardisoSolver(const SparseMatrix * A, int numThreads, int positiveDefinite=0, int directIterative=0, int verbose=0);
@@ -87,7 +87,7 @@ protected:
void *pt[64];
MKL_INT iparm[64];
int mtype;
- MKL_INT nrhs;
+ MKL_INT nrhs;
MKL_INT maxfct, mnum, phase, error, msglvl;
static void DisabledSolverError();
diff --git a/src/libsparseSolver/SPOOLESSolver.cpp b/src/libsparseSolver/SPOOLESSolver.cpp
index 24ee33349585e869d9ec28324672f0dd0432150e..812ef9a0ee1ba8d5018456358222b205c8da5e34 100644
--- a/src/libsparseSolver/SPOOLESSolver.cpp
+++ b/src/libsparseSolver/SPOOLESSolver.cpp
@@ -36,9 +36,9 @@
unlike, say, in the Conjugate gradient method.
Memory requirements are minimized by re-ordering the matrix before applying
Cholesky decomposition.
- However, for very large systems (e.g. 200,000 x 200,000 matrices on a
+ However, for very large systems (e.g. 200,000 x 200,000 matrices on a
2Gb machine), the Cholesky decomposition might run out of memory.
-
+
Jernej Barbic, MIT, 2007-2009
*/
@@ -152,7 +152,7 @@ int SPOOLESSolver::SolveLinearSystem(double * x, const double * rhs)
// zero the solution
DenseMtx_zero(mtx_x);
-
+
// solve
if (verbose >= 2)
printf("Solving the linear system...\n");
@@ -167,7 +167,7 @@ int SPOOLESSolver::SolveLinearSystem(double * x, const double * rhs)
}
if (verbose >= 2)
- printf("Solve completed.\n");
+ printf("Solve completed.\n");
// store result
for(int i=0; i < n; i++)
diff --git a/src/libsparseSolver/SPOOLESSolver.h b/src/libsparseSolver/SPOOLESSolver.h
index 40f229adc8c673c4786f9d5a2a3c0e2e3f827f0c..124f279b5f7d5cd204cfe16989eb90d8dfeadd37 100644
--- a/src/libsparseSolver/SPOOLESSolver.h
+++ b/src/libsparseSolver/SPOOLESSolver.h
@@ -32,16 +32,16 @@
/*
Solves A * x = rhs, where A is sparse, usually large, and symmetric.
-
+
The solution is obtained using the SPOOLES library (which is free software).
The solution method is direct (not iterative). As such, convergence
is often very robust, and there is no need to tune convergence parameters,
unlike, say, in the Conjugate gradient method.
Memory requirements are minimized by re-ordering the matrix before applying
Cholesky decomposition.
- However, for very large systems (e.g. 200,000 x 200,000 matrices on a
+ However, for very large systems (e.g. 200,000 x 200,000 matrices on a
2Gb machine), the Cholesky decomposition might run out of memory.
-
+
Jernej Barbic, MIT, 2007-2009
*/
diff --git a/src/libsparseSolver/SPOOLESSolverMT.h b/src/libsparseSolver/SPOOLESSolverMT.h
index 410edbe923ad578640fea6d987b90fe9895f3306..d88befb3f2a1d6ebde60ee2ed694381eb19e9fcc 100644
--- a/src/libsparseSolver/SPOOLESSolverMT.h
+++ b/src/libsparseSolver/SPOOLESSolverMT.h
@@ -40,7 +40,7 @@ is often very robust, and there is no need to tune convergence parameters,
unlike, say, in the Conjugate gradient method.
Memory requirements are minimized by re-ordering the matrix before applying
Cholesky decomposition.
-However, for very large systems (e.g. 200,000 x 200,000 matrices on a
+However, for very large systems (e.g. 200,000 x 200,000 matrices on a
2Gb machine), the Cholesky decomposition might run out of memory.
*/
diff --git a/src/libsparseSolver/linearSolver.h b/src/libsparseSolver/linearSolver.h
index bd2430e981baecec2e50b72cedb3125343576533..f17eb200e0a0fe2dad40029fc9880437092279bf 100644
--- a/src/libsparseSolver/linearSolver.h
+++ b/src/libsparseSolver/linearSolver.h
@@ -32,7 +32,7 @@
/*
Abstract class to solve
A * x = rhs, where A is a square matrix.
-
+
Jernej Barbic, USC, 2010
*/
diff --git a/src/util/interactiveDeformableSimulator/interactiveDeformableSimulator.cpp b/src/util/interactiveDeformableSimulator/interactiveDeformableSimulator.cpp
index a2cd566f30924e24f3f38c4ec73758979663202a..c3176f71de6772d7b811c73104e57c98a30265a9 100755
--- a/src/util/interactiveDeformableSimulator/interactiveDeformableSimulator.cpp
+++ b/src/util/interactiveDeformableSimulator/interactiveDeformableSimulator.cpp
@@ -1509,7 +1509,6 @@ void initSimulation()
{
case INV_STVK:
{
-
isotropicMaterial = new StVKIsotropicMaterial(tetMesh, enableCompressionResistance, compressionResistance);
printf("Invertible material: StVK.\n");
break;