*/
class dSMatrix
{
- dSMatrixStruct& cols;
- double* val;
- int max_len;
- public:
- //////////
- void reinit();
- //////////
- dSMatrix(dSMatrixStruct& c)
- : cols(c), val(0), max_len(0)
- {
- reinit();
- }
- //////////
- ~dSMatrix()
- {
- delete[] val;
- }
+ dSMatrixStruct * cols;
+ double* val;
+ int max_len;
+ public:
+ //
+ void reinit();
+ //
+ void reinit (dSMatrixStruct &);
+ //
+ dSMatrix(dSMatrixStruct& c)
+ : cols(&c), val(0), max_len(0)
+ {
+ reinit();
+ }
+ //
+ ~dSMatrix()
+ {
+ delete[] val;
+ }
- //////////
- int m() const { return cols.rows; }
- //////////
- int n() const { return cols.cols; }
+ //
+ int m() const { return cols->rows; }
+ //
+ int n() const { return cols->cols; }
- //////////
- void set(int i,int j,double value) { val[cols(i,j)] = value; }
- //////////
- void add(int i,int j,double value) { val[cols(i,j)]+= value; }
+ //
+ void set(int i,int j,double value) { val[cols->operator()(i,j)] = value; }
+ //
+ void add(int i,int j,double value) { val[cols->operator()(i,j)]+= value; }
- //////////
- void vmult (dVector& dst,const dVector& src);
- //////////
- void Tvmult(dVector& dst,const dVector& src);
+ //
+ void vmult (dVector& dst,const dVector& src);
+ //
+ void Tvmult(dVector& dst,const dVector& src);
- //////////
- double residual (dVector& dst,const dVector& x,const dVector& b);
- //////////
- void Jacobi_precond(dVector& dst,const dVector& src,double om = 1.);
- //////////
- void SSOR_precond(dVector& dst,const dVector& src,double om = 1.);
- //////////
- void SOR_precond(dVector& dst,const dVector& src,double om = 1.);
- //////////
- void SSOR(dVector& dst,double om = 1.);
- //////////
- void SOR(dVector& dst,double om = 1.);
- //////////
- void precondition(dVector& dst,const dVector& src) { dst=src; }
+ //
+ double residual (dVector& dst,const dVector& x,const dVector& b);
+ //
+ void Jacobi_precond(dVector& dst,const dVector& src,double om = 1.);
+ //
+ void SSOR_precond(dVector& dst,const dVector& src,double om = 1.);
+ //
+ void SOR_precond(dVector& dst,const dVector& src,double om = 1.);
+ //
+ void SSOR(dVector& dst,double om = 1.);
+ //
+ void SOR(dVector& dst,double om = 1.);
+ //
+ void precondition(dVector& dst,const dVector& src) { dst=src; }
void print (ostream &) const;
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcNotCompressed);
};
#endif
add(rows(i), cols(j));
}
-void
-dSMatrix::reinit()
-{
- if(max_len<cols.vec_len)
- {
- if (val) delete[] val;
- val = new double[cols.vec_len];
- max_len = cols.vec_len;
- }
-// memset(val, 0, sizeof(*val) * cols.vec_len);
- for (int i = cols.vec_len-1 ; i>=0 ; i--) val[i] = 0;
-}
void
dSMatrixStruct::print_gnuplot (ostream &out) const
return b;
}
+
+
+void
+dSMatrix::reinit()
+{
+ Assert (cols->compressed, ExcNotCompressed());
+
+ if(max_len<cols->vec_len)
+ {
+ if (val) delete[] val;
+ val = new double[cols->vec_len];
+ max_len = cols->vec_len;
+ }
+// memset(val, 0, sizeof(*val) * cols->vec_len);
+ for (int i = cols->vec_len-1 ; i>=0 ; i--) val[i] = 0;
+}
+
+
+
+void
+dSMatrix::reinit (dSMatrixStruct &sparsity) {
+ cols = &sparsity;
+ reinit ();
+};
+
+
+
+
void
dSMatrix::vmult(dVector& dst,const dVector& src)
{
for (int i=0;i<m();i++)
{
double s = 0.;
- for (int j=cols.rowstart[i]; j < cols.rowstart[i+1] ;j++)
+ for (int j=cols->rowstart[i]; j < cols->rowstart[i+1] ;j++)
{
- int p = cols.colnums[j];
+ int p = cols->colnums[j];
s += val[j] * src(p);
}
dst(i) = s;
for (i=0;i<m();i++)
{
- for (int j=cols.rowstart[i]; j<cols.rowstart[i+1] ;j++)
+ for (int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
- int p = cols.colnums[j];
+ int p = cols->colnums[j];
dst(p) += val[j] * src(i);
}
}
for (int i=0;i<m();i++)
{
s = b(i);
- for (int j=cols.rowstart[i]; j<cols.rowstart[i+1] ;j++)
+ for (int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
- int p = cols.colnums[j];
+ int p = cols->colnums[j];
s -= val[j] * u(p);
}
dst(i) = s;
for (int i=0;i<n;++i)
{
- dst(i) = om * src(i) * val[cols.rowstart[i]];
+ dst(i) = om * src(i) * val[cols->rowstart[i]];
}
}
for (i=0;i<n;i++)
{
dst(i) = src(i);
- for (j=cols.rowstart[i]; j<cols.rowstart[i+1] ;j++)
+ for (j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
- p = cols.colnums[j];
+ p = cols->colnums[j];
if (p<i) dst(i) -= om* val[j] * dst(p);
}
- dst(i) /= val[cols.rowstart[i]];
+ dst(i) /= val[cols->rowstart[i]];
}
- for (i=0;i<n;i++) dst(i) *= (2.-om)*val[cols.rowstart[i]];
+ for (i=0;i<n;i++) dst(i) *= (2.-om)*val[cols->rowstart[i]];
for (i=n-1;i>=0;i--)
{
- for (j=cols.rowstart[i];j<cols.rowstart[i+1];j++)
+ for (j=cols->rowstart[i];j<cols->rowstart[i+1];j++)
{
- p = cols.colnums[j];
+ p = cols->colnums[j];
if (p>i) dst(i) -= om* val[j] * dst(p);
}
- dst(i) /= val[cols.rowstart[i]];
+ dst(i) /= val[cols->rowstart[i]];
}
}
for (int i=0;i<m();i++)
{
double s = dst(i);
- for (int j=cols.rowstart[i]; j<cols.rowstart[i+1] ;j++)
+ for (int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
- int p = cols.colnums[j];
+ int p = cols->colnums[j];
if (p<i)
s -= val[j] * dst(p);
}
- dst(i) = s * om / val[cols.rowstart[i]];
+ dst(i) = s * om / val[cols->rowstart[i]];
}
}
for (i=0;i<n;i++)
{
s = 0.;
- for (j=cols.rowstart[i]; j<cols.rowstart[i+1] ;j++)
+ for (j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
- p = cols.colnums[j];
+ p = cols->colnums[j];
if (p>=0)
{
if (i>j) s += val[j] * dst(p);
}
}
dst(i) -= s * om;
- dst(i) /= val[cols.rowstart[i]];
+ dst(i) /= val[cols->rowstart[i]];
}
for (i=n-1;i>=0;i--)
{
s = 0.;
- for (j=cols.rowstart[i]; j<cols.rowstart[i+1] ;j++)
+ for (j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
- p = cols.colnums[j];
+ p = cols->colnums[j];
if (p>=0)
{
if (i<j) s += val[j] * dst(p);
}
}
- dst(i) -= s * om / val[cols.rowstart[i]];
+ dst(i) -= s * om / val[cols->rowstart[i]];
}
}
void dSMatrix::print (ostream &out) const {
- for (int i=0; i<cols.rows; ++i)
- for (int j=cols.rowstart[i]; j<cols.rowstart[i+1]; ++j)
- out << "(" << i << "," << cols.colnums[j] << ") " << val[j] << endl;
+ for (int i=0; i<cols->rows; ++i)
+ for (int j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
+ out << "(" << i << "," << cols->colnums[j] << ") " << val[j] << endl;
};