dolfin team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

On 8 Mar 2012, at 21:43, Johan Hake <hake.dev@xxxxxxxxx> wrote:

> On Thu, Mar 8, 2012 at 9:40 PM, Garth N. Wells <gnw20@xxxxxxxxx> wrote:
>> This change is not safe. For a matrix and a vector to be compatible in
>> parallel, it is not sufficient to check the global size.
> 
> True. The logic is that a user need to provide a result vector with
> correct size and parallell structure. If not just pass an empty
> vector. The problem with always resizing the result is that one loose
> performance for iterative solvers using A.mult(x, y).


The only proper solution that I see is that a user is responsible for making sure that the size is right and that the function do no resizing. The change is a half-way house, so I think that it should be reverted. The old code was safe, and the proper, efficient alternative is that no resizing be done.

Garth

> 
> Johan
> 
>> Garth
>> 
>> 
>> 
>> 
>> ---------- Forwarded message ----------
>> From:  <noreply@xxxxxxxxxxxxx>
>> Date: 8 March 2012 20:31
>> Subject: [Branch ~dolfin-core/dolfin/trunk] Rev 6630: Avoid
>> unnessesary resize of result vector for A.mult(x, y)
>> To: Garth Wells <gnw20@xxxxxxxxx>
>> 
>> 
>> ------------------------------------------------------------
>> revno: 6630
>> committer: Johan Hake <hake.dev@xxxxxxxxx>
>> branch nick: work-trunk
>> timestamp: Thu 2012-03-08 21:25:41 +0100
>> message:
>>  Avoid unnessesary resize of result vector for A.mult(x, y)
>>    -- Apply patch from Joachim
>> modified:
>>  ChangeLog
>>  dolfin/la/EpetraMatrix.cpp
>>  dolfin/la/GenericMatrix.h
>>  dolfin/la/PETScMatrix.cpp
>>  dolfin/la/uBLASMatrix.h
>>  dolfin/swig/la/post.i
>>  test/unit/la/python/KrylovSolver.py
>> 
>> 
>> --
>> lp:dolfin
>> https://code.launchpad.net/~dolfin-core/dolfin/trunk
>> 
>> Your team DOLFIN Core Team is subscribed to branch lp:dolfin.
>> To unsubscribe from this branch go to
>> https://code.launchpad.net/~dolfin-core/dolfin/trunk/+edit-subscription
>> 
>> === modified file 'ChangeLog'
>> --- ChangeLog   2012-02-13 21:49:10 +0000
>> +++ ChangeLog   2012-03-08 20:25:41 +0000
>> @@ -1,3 +1,4 @@
>> + - Avoid unnessesary resize of result vector for A*b
>>  - Major speed-up of mesh topology computation
>>  - Add simple 2D and 3D mesh generation (via CGAL)
>>  - Add creation of mesh from triangulations of points (via CGAL)
>> 
>> === modified file 'dolfin/la/EpetraMatrix.cpp'
>> --- dolfin/la/EpetraMatrix.cpp  2012-02-24 12:24:53 +0000
>> +++ dolfin/la/EpetraMatrix.cpp  2012-03-08 20:25:41 +0000
>> @@ -568,7 +568,17 @@
>>   }
>> 
>>   // Resize RHS
>> -  this->resize(Ax, 0);
>> +  if (Ax.size() == 0)
>> +  {
>> +    this->resize(Ax, 0);
>> +  }
>> +
>> +  if (Ax.size() != size(0))
>> +  {
>> +    dolfin_error("EpetraMatrix.cpp",
>> +                 "compute matrix-vector product with Epetra matrix",
>> +                 "Vector for matrix-vector result has wrong size");
>> +  }
>> 
>>   dolfin_assert(x.vec());
>>   dolfin_assert(Ax.vec());
>> @@ -595,7 +605,17 @@
>>   }
>> 
>>   // Resize RHS
>> -  this->resize(Ax, 1);
>> +  if (Ax.size() == 0)
>> +  {
>> +    this->resize(Ax, 1);
>> +  }
>> +
>> +  if (Ax.size() != size(1))
>> +  {
>> +    dolfin_error("EpetraMatrix.cpp",
>> +                 "compute transpose matrix-vector product with Epetra matrix",
>> +                 "Vector for transpose matrix-vector result has wrong size");
>> +  }
>> 
>>   const int err = A->Multiply(true, *(x.vec()), *(Ax.vec()));
>>   if (err != 0)
>> 
>> === modified file 'dolfin/la/GenericMatrix.h'
>> --- dolfin/la/GenericMatrix.h   2012-02-28 22:23:00 +0000
>> +++ dolfin/la/GenericMatrix.h   2012-03-08 20:25:41 +0000
>> @@ -137,10 +137,12 @@
>>     /// Set given rows to identity matrix
>>     virtual void ident(uint m, const uint* rows) = 0;
>> 
>> -    /// Matrix-vector product, y = Ax
>> +    /// Matrix-vector product, y = Ax. The y vector must either be zero-sized
>> +    /// or have correct size and parallel layout.
>>     virtual void mult(const GenericVector& x, GenericVector& y) const = 0;
>> 
>> -    /// Matrix-vector product, y = A^T x
>> +    /// Matrix-vector product, y = A^T x. The y vector must either be
>> +    /// zero-sized or have correct size and parallel layout.
>>     virtual void transpmult(const GenericVector& x, GenericVector& y) const = 0;
>> 
>>     /// Multiply matrix by given number
>> 
>> === modified file 'dolfin/la/PETScMatrix.cpp'
>> --- dolfin/la/PETScMatrix.cpp   2012-02-24 12:24:53 +0000
>> +++ dolfin/la/PETScMatrix.cpp   2012-03-08 20:25:41 +0000
>> @@ -393,8 +393,19 @@
>>                  "Non-matching dimensions for matrix-vector product");
>>   }
>> 
>> -  // Resize if dimensions does not match
>> -  resize(yy, 0);
>> +  // Resize RHS if empty
>> +  if (yy.size() == 0)
>> +  {
>> +    resize(yy, 0);
>> +  }
>> +
>> +  if (size(0) != yy.size())
>> +  {
>> +    dolfin_error("PETScMatrix.cpp",
>> +                 "compute matrix-vector product with PETSc matrix",
>> +                 "Vector for matrix-vector result has wrong size");
>> +  }
>> +
>>   MatMult(*A, *xx.vec(), *yy.vec());
>>  }
>>  //-----------------------------------------------------------------------------
>> @@ -412,7 +423,19 @@
>>                  "Non-matching dimensions for transpose matrix-vector product");
>>   }
>> 
>> -  resize(yy, 1);
>> +  // Resize RHS if empty
>> +  if (yy.size() == 0)
>> +  {
>> +    resize(yy, 1);
>> +  }
>> +
>> +  if (size(1) != yy.size())
>> +  {
>> +    dolfin_error("PETScMatrix.cpp",
>> +                 "compute transpose matrix-vector product with PETSc matrix",
>> +                 "Vector for transpose matrix-vector result has wrong size");
>> +  }
>> +
>>   MatMultTranspose(*A, *xx.vec(), *yy.vec());
>>  }
>>  //-----------------------------------------------------------------------------
>> 
>> === modified file 'dolfin/la/uBLASMatrix.h'
>> --- dolfin/la/uBLASMatrix.h     2012-02-24 12:24:53 +0000
>> +++ dolfin/la/uBLASMatrix.h     2012-03-08 20:25:41 +0000
>> @@ -461,7 +461,31 @@
>>   template <typename Mat>
>>   void uBLASMatrix<Mat>::mult(const GenericVector& x, GenericVector& y) const
>>   {
>> -    ublas::axpy_prod(A, x.down_cast<uBLASVector>().vec(),
>> y.down_cast<uBLASVector>().vec(), true);
>> +
>> +    const uBLASVector& xx = x.down_cast<uBLASVector>();
>> +    uBLASVector& yy = y.down_cast<uBLASVector>();
>> +
>> +    if (size(1) != xx.size())
>> +    {
>> +      dolfin_error("uBLASMatrix.cpp",
>> +                   "compute matrix-vector product with uBLAS matrix",
>> +                   "Non-matching dimensions for matrix-vector product");
>> +    }
>> +
>> +  // Resize RHS if empty
>> +  if (yy.size() == 0)
>> +  {
>> +    resize(yy, 0);
>> +  }
>> +
>> +  if (size(0) != yy.size())
>> +  {
>> +    dolfin_error("uBLASMatrix.cpp",
>> +                 "compute matrix-vector product with uBLAS matrix",
>> +                 "Vector for matrix-vector result has wrong size");
>> +  }
>> +
>> +    ublas::axpy_prod(A, xx.vec(), yy.vec(), true);
>>   }
>>   //-----------------------------------------------------------------------------
>>   template <typename Mat>
>> 
>> === modified file 'dolfin/swig/la/post.i'
>> --- dolfin/swig/la/post.i       2012-02-21 21:45:42 +0000
>> +++ dolfin/swig/la/post.i       2012-03-08 20:25:41 +0000
>> @@ -556,9 +556,9 @@
>>             if vector_type not in _matrix_vector_mul_map[matrix_type]:
>>                 raise TypeError, "Provide a Vector which can be
>> down_casted to ''"%vector_type.__name__
>>             if type(other) == Vector:
>> -                ret = Vector(self.size(0))
>> +                ret = Vector()
>>             else:
>> -                ret = vector_type(self.size(0))
>> +                ret = vector_type()
>>             self.mult(other, ret)
>>             return ret
>>         elif isinstance(other, ndarray):
>> 
>> === modified file 'test/unit/la/python/KrylovSolver.py'
>> --- test/unit/la/python/KrylovSolver.py 2012-02-29 18:44:36 +0000
>> +++ test/unit/la/python/KrylovSolver.py 2012-03-08 20:25:41 +0000
>> @@ -84,15 +84,12 @@
>>                 def __init__(self, A):
>>                     self.A = A
>> 
>> -                    # Use tmp vector as A.mult will resize the vector
>> -                    self.y_tmp = PETScVector()
>>                     PETScKrylovMatrix.__init__(self, A.size(0), A.size(1))
>> 
>>                 def mult(self, x, y):
>> 
>>                     # Make ordinary matrix vector product
>> -                    self.A.mult(x, self.y_tmp)
>> -                    y[:] = self.y_tmp
>> +                    self.A.mult(x, y)
>> 
>> 
>>             class MatrixFreeKrylovMatrix(PETScKrylovMatrix) :
>> 
>> _______________________________________________
>> Mailing list: https://launchpad.net/~dolfin
>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>> Unsubscribe : https://launchpad.net/~dolfin
>> More help   : https://help.launchpad.net/ListHelp
>>