solve Px = rhs
*/
Vector solve(Vector rhs) const;
-
+ void solve (Vector &dest, Vector const &rhs)const;
Vector operator * (Vector rhs) const { return solve (rhs); }
/**
return the inverse of the matrix P.
*/
Matrix original() const;
private:
+ void full_matrix_solve(Vector &,Vector const&)const;
+ void band_matrix_solve(Vector &, Vector const&)const;
void full_matrix_decompose(Matrix const & P);
void band_matrix_decompose(Matrix const &P);
int
Matrix::band_i()const
{
+ if ( band_b() ) {
+ Diagonal_storage const * diag = (Diagonal_storage*) dat;
+ return diag->band_size_i();
+ }
int starty = dim();
while (starty >= 0 ) {
for ( int i = starty, j = 0; i < dim(); i++, j++ )
return m;
}
+/*
+ ugh. Only works for square matrices.
+ */
void
Matrix::transpose() // delegate to storage?
{
+#if 1
for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) {
if (i >= j)
continue;
dat->elem(i,j) = dat->elem(j,i);
dat->elem(j,i)=r;
}
+#endif
}
Matrix
if (m1.dat->try_right_multiply(dat, m2.dat))
return;
+
for (int i=0, j=0; dat->mult_ok(i,j);
dat->mult_next(i,j)) {
Real r=0.0;