- petsc-dev - mailman3.cels.anl.gov

Re: patch for seqdense
by S V N Vishwanathan 28 Aug '05

28 Aug '05

Hi! >> I sent through a patch to implement MatMatMult* for seqdense matrices. Hong> Where is this patch? through bk or an attached file? I don't see Hong> MatMatMult_seqdense from Hong> petsc-dev. Hong> Also, do you have testing code for it? I can test it Hong> once I see the MatMatMult_seqdense. I think Simon forgot to add the patch file. Here is it attached again. Please let me know if it applies cleanly. We do not have bk but can consider licensing it if PETSc is committed to using it in the near future and we can get write access to the source code. But, we would much rather prefer to use an open source SCM tool. vishy diff -ur petsc-bk/src/mat/impls/dense/mpi/makefile petsc-dev/src/mat/impls/dense/mpi/makefile --- petsc-bk/src/mat/impls/dense/mpi/makefile 2005-08-17 10:34:32.000000000 +1000 +++ petsc-dev/src/mat/impls/dense/mpi/makefile 2005-08-17 10:38:12.000000000 +1000 @@ -9,7 +9,7 @@ OBJSC = mpidense.o mmdense.o OBJSF = LIBBASE = libpetscmat -DIRS = +DIRS = plapack MANSEC = Mat LOCDIR = src/mat/impls/dense/mpi/ diff -ur petsc-bk/src/mat/impls/dense/seq/dense.c petsc-dev/src/mat/impls/dense/seq/dense.c --- petsc-bk/src/mat/impls/dense/seq/dense.c 2005-08-16 14:41:57.000000000 +1000 +++ petsc-dev/src/mat/impls/dense/seq/dense.c 2005-08-17 13:18:12.000000000 +1000 @@ -1319,7 +1319,7 @@ Mat_SeqDense *mat = (Mat_SeqDense*)A->data; PetscFunctionBegin; - if (mat->lda != A->m) SETERRQ(PETSC_ERR_SUP,"Cannot get array for Denses matrices with LDA different from number of rows"); + if (mat->lda != A->m) SETERRQ(PETSC_ERR_SUP,"Cannot get array for Dense matrices with LDA different from number of rows"); *array = mat->v; PetscFunctionReturn(0); } @@ -1476,6 +1476,94 @@ PetscFunctionReturn(0); } +/* ----------------------------------------------------------------*/ + +#undef __FUNCT__ +#define __FUNCT__ "MatMatMult_SeqDense_SeqDense" +PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) +{ + PetscErrorCode ierr; + + PetscFunctionBegin; + if (scall == MAT_INITIAL_MATRIX){ + ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); + } + ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); + PetscFunctionReturn(0); +} + + +#undef __FUNCT__ +#define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" +PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) +{ + PetscFunctionBegin; + PetscFunctionReturn(PETSC_ERR_SUP); +} + + +#undef __FUNCT__ +#define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" +PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) +{ + Mat_SeqDense *a = (Mat_SeqDense*)A->data; + Mat_SeqDense *b = (Mat_SeqDense*)B->data; + Mat_SeqDense *c = (Mat_SeqDense*)C->data; + PetscBLASInt m = (PetscBLASInt)A->m, n = (PetscBLASInt)B->n, k = (PetscBLASInt)A->n; + PetscScalar _DOne=1.0,_DZero=0.0; + /* PetscErrorCode ierr; */ + + PetscFunctionBegin; + + BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); + + PetscFunctionReturn(0); +} + + +#undef __FUNCT__ +#define __FUNCT__ "MatMatMultTranspose_SeqDense_SeqDense" +PetscErrorCode MatMatMultTranspose_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) +{ + PetscErrorCode ierr; + + PetscFunctionBegin; + if (scall == MAT_INITIAL_MATRIX){ + ierr = MatMatMultTransposeSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); + } + ierr = MatMatMultTransposeNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); + PetscFunctionReturn(0); +} + + +#undef __FUNCT__ +#define __FUNCT__ "MatMatMultTransposeSymbolic_SeqDense_SeqDense" +PetscErrorCode MatMatMultTransposeSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) +{ + PetscFunctionBegin; + PetscFunctionReturn(PETSC_ERR_SUP); +} + + +#undef __FUNCT__ +#define __FUNCT__ "MatMatMultTransposeNumeric_SeqDense_SeqDense" +PetscErrorCode MatMatMultTransposeNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) +{ + Mat_SeqDense *a = (Mat_SeqDense*)A->data; + Mat_SeqDense *b = (Mat_SeqDense*)B->data; + Mat_SeqDense *c = (Mat_SeqDense*)C->data; + PetscBLASInt m = (PetscBLASInt)A->n, n = (PetscBLASInt)B->n, k = (PetscBLASInt)A->m; + PetscScalar _DOne=1.0,_DZero=0.0; + /* PetscErrorCode ierr; */ + + PetscFunctionBegin; + + BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); + + PetscFunctionReturn(0); +} + + /* -------------------------------------------------------------------*/ static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, MatGetRow_SeqDense, @@ -1566,16 +1654,16 @@ 0, 0, 0, - 0, -/*90*/ 0, - 0, 0, +/*90*/ MatMatMult_SeqDense_SeqDense, + MatMatMultSymbolic_SeqDense_SeqDense, + MatMatMultNumeric_SeqDense_SeqDense, 0, 0, /*95*/ 0, - 0, - 0, - 0, + MatMatMultTranspose_SeqDense_SeqDense, + MatMatMultTransposeSymbolic_SeqDense_SeqDense, + MatMatMultTransposeNumeric_SeqDense_SeqDense, 0, /*100*/0, 0, diff -ur petsc-bk/src/mat/impls/dense/seq/dense.h petsc-dev/src/mat/impls/dense/seq/dense.h --- petsc-bk/src/mat/impls/dense/seq/dense.h 2005-07-11 12:46:24.000000000 +1000 +++ petsc-dev/src/mat/impls/dense/seq/dense.h 2005-08-17 11:21:09.000000000 +1000 @@ -23,5 +23,11 @@ EXTERN PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec,Vec,Vec); EXTERN PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec,Vec); EXTERN PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec,Vec,Vec); +EXTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat,Mat,MatReuse,PetscReal,Mat*); +EXTERN PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat,Mat,PetscReal,Mat*); +EXTERN PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat,Mat,Mat); +EXTERN PetscErrorCode MatMatMultTranspose_SeqDense_SeqDense(Mat,Mat,MatReuse,PetscReal,Mat*); +EXTERN PetscErrorCode MatMatMultTransposeSymbolic_SeqDense_SeqDense(Mat,Mat,PetscReal,Mat*); +EXTERN PetscErrorCode MatMatMultTransposeNumeric_SeqDense_SeqDense(Mat,Mat,Mat); #endif ----- But the possibility of abuse may be a good reason for leaving capabilities out of other computer languages, it's not a good reason for leaving capabilities out of Perl. -- Larry Wall in <199709251614.JAA15718(a)wall.org>

1 0

Re: patch for seqdense
by Hong Zhang 28 Aug '05

28 Aug '05

Simon, > > I sent through a patch to implement MatMatMult* for seqdense matrices. Where is this patch? through bk or an attached file? I don't see MatMatMult_seqdense from petsc-dev. Also, do you have testing code for it? I can test it once I see the MatMatMult_seqdense. Thanks, Hong > Here it is again; is there anything wrong with it ? It just calls through to BLASgemm_ . > > thanks, > > Simon. > > -- > Simon Burton, B.Sc. > Licensed PO Box 8066 > ANU Canberra 2601 > Australia > Ph. 61 02 6249 6940 > http://arrowtheory.com > >

1 0

patch for seqdense
by Simon Burton 27 Aug '05

27 Aug '05

Hi, I sent through a patch to implement MatMatMult* for seqdense matrices. Here it is again; is there anything wrong with it ? It just calls through to BLASgemm_ . thanks, Simon. -- Simon Burton, B.Sc. Licensed PO Box 8066 ANU Canberra 2601 Australia Ph. 61 02 6249 6940 http://arrowtheory.com

1 0

Re: Fortran interfaces
by William Gropp 22 Aug '05

22 Aug '05

At 09:06 AM 8/19/2005, Barry Smith wrote: > Bill, > > Do you think it would be possible to enhance bfort >to generate Fortran 90 interface definitions? Something like > > PetscErrorCode MatFactorInfoInitialize(MatFactorInfo info) > >becomes > > interface > subroutine MatFactorInfoInitialize(info,ierr) > MatFactorInfo :: info(11) > PetscErrorCode :: ierr > end subroutine > end interface > > I think this would help with Fortran programmers writing correct >codes much quicker. I think that this is a good idea. This would be a good student project. Bill > Thanks, > > Barry William Gropp http://www.mcs.anl.gov/~gropp

1 0

Re: ufuncs, iterators
by Barry Smith 20 Aug '05

20 Aug '05

I withdraw my concern for the (I thought) possible mislocation of the various operations on "Matrices". So folks should feel free to begin to add them. Barry On Sat, 20 Aug 2005, S V N Vishwanathan wrote: > Hi! > > Barry> After sleeping on it, it may be ok to include these methods in > Barry> Mat as methods for "filling up matrices" so long as in the end > Barry> you end up with a Mat that you then use as an operator. > > Yes that sounds like a good test to me. If the end result of the > computation is a matrix which is going to be used as a linear operator > then we must include it in. That means that for instance methods like > exp(M) and M.*X should also be supported. Right? > > Barry> But I'd still like to see/understand a little more of the > Barry> "construction" process. Classically one would do that as loops > Barry> over elements and perform all the computations for the one > Barry> element before moving to the next. In Matlab this can be done > Barry> instead (with some impact on performance) using a sequence of > Barry> array operations (with the loop inside each array operation). In > Barry> the past, since PETSc was used exclusively from C/C++ and Fortran > Barry> users built their matrices (operators) using the "classical" > Barry> approach, now with python it appears reasonable that we may need > Barry> to add the "array" approach. > > Yes I think we should add both i.e. ability to loop over all elements > and apply a function or to loop only over the arrays and perform > operations on the arrays. > > >> I am thinking more about what Barry said. The VecPointwise*() > >> operations can be given a solid mathematical interpretation in terms > >> of spinor operations. However, I do not see anything like that for > >> the Mat stuff yet. We need to understand the mathematicas better. > > Pardon my ignorance. What are spinor operations? Why isn't something > like a M.*X a matrix operation. It is the Hadamard product. > > vishy > > ----- > Don't you wish that all the people who sincerely want to help you > could agree with each other? > >

1 0

Re: ufuncs, iterators
by Simon Burton 20 Aug '05

20 Aug '05

Oh, that should probably read: exp( -1/2\sigma^{2} ||x1_{i} - x2_{j}||_{2}^{2}) And when we vectorize this operation: ||x1_i - x2_j||^2 = ||x1_i||^2 + ||x2_i||^2 + 2*(x1_i,x2_j) and the last term is the ip matrix. It seems that this expression satisfies the linear operator requirement (it's invarient under linear isometries). Simon. On Thu, 18 Aug 2005 23:12:21 -0500 (CDT) Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > In terms of exp( -1/2\sigma^{2} ||x_{i} - x_{j}||_{2}^{2}) > what are they? > > Thanks > > Barry > > > On Fri, 19 Aug 2005, Simon Burton wrote: > > > On Thu, 18 Aug 2005 22:49:29 -0500 (CDT) > > Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > > > > > > What is x1, x2 and ip? > > > > > > Barry > > > > x1 and x2 are 2-arrays; their rows are the 'sample' vectors. > > ip is the matrix of all inner products from x1 and x2. > > > > Simon. > > > > > > > -- Simon Burton, B.Sc. Licensed PO Box 8066 ANU Canberra 2601 Australia Ph. 61 02 6249 6940 http://arrowtheory.com

1 0

Re: ufuncs, iterators
by Matthew Knepley 20 Aug '05

20 Aug '05

Simon Burton <simon(a)arrowtheory.com> writes: > Oh, that should probably read: > > exp( -1/2\sigma^{2} ||x1_{i} - x2_{j}||_{2}^{2}) > > And when we vectorize this operation: > > ||x1_i - x2_j||^2 = ||x1_i||^2 + ||x2_i||^2 + 2*(x1_i,x2_j) > > and the last term is the ip matrix. > > It seems that this expression satisfies the linear operator > requirement (it's invarient under linear isometries). I am looking again, and it appears to be the normalized graph laplacian. This is indeed a linear operator which we need to support. Matt > Simon. > > > On Thu, 18 Aug 2005 23:12:21 -0500 (CDT) > Barry Smith <bsmith(a)mcs.anl.gov> wrote: > >> >> In terms of exp( -1/2\sigma^{2} ||x_{i} - x_{j}||_{2}^{2}) >> what are they? >> >> Thanks >> >> Barry >> >> >> On Fri, 19 Aug 2005, Simon Burton wrote: >> >> > On Thu, 18 Aug 2005 22:49:29 -0500 (CDT) >> > Barry Smith <bsmith(a)mcs.anl.gov> wrote: >> > >> > > >> > > What is x1, x2 and ip? >> > > >> > > Barry >> > >> > x1 and x2 are 2-arrays; their rows are the 'sample' vectors. >> > ip is the matrix of all inner products from x1 and x2. >> > >> > Simon. >> > >> > >> > >> > > > -- > Simon Burton, B.Sc. > Licensed PO Box 8066 > ANU Canberra 2601 > Australia > Ph. 61 02 6249 6940 > http://arrowtheory.com > > > -- "Failure has a thousand explanations. Success doesn't need one" -- Sir Alec Guiness

1 0

Re: ufuncs, iterators
by S V N Vishwanathan 20 Aug '05

20 Aug '05

Hi! Barry> Certainly pointwise addition is the same as adding two operators Barry> together (PETSc has this with MatAXPY), but what about, for example, Barry> pointwise multiply? Pointwise multiply is the Hadamard product of two matrices matrix. You have something like MatGetRowMax which strictly speaking is not a matrix operation in the linear algebra sense. So why not something like MatGetRowNorm? vishy

1 0

Re: ufuncs, iterators
by S V N Vishwanathan 20 Aug '05

20 Aug '05

Hi! Barry> After sleeping on it, it may be ok to include these methods in Barry> Mat as methods for "filling up matrices" so long as in the end Barry> you end up with a Mat that you then use as an operator. Yes that sounds like a good test to me. If the end result of the computation is a matrix which is going to be used as a linear operator then we must include it in. That means that for instance methods like exp(M) and M.*X should also be supported. Right? Barry> But I'd still like to see/understand a little more of the Barry> "construction" process. Classically one would do that as loops Barry> over elements and perform all the computations for the one Barry> element before moving to the next. In Matlab this can be done Barry> instead (with some impact on performance) using a sequence of Barry> array operations (with the loop inside each array operation). In Barry> the past, since PETSc was used exclusively from C/C++ and Fortran Barry> users built their matrices (operators) using the "classical" Barry> approach, now with python it appears reasonable that we may need Barry> to add the "array" approach. Yes I think we should add both i.e. ability to loop over all elements and apply a function or to loop only over the arrays and perform operations on the arrays. >> I am thinking more about what Barry said. The VecPointwise*() >> operations can be given a solid mathematical interpretation in terms >> of spinor operations. However, I do not see anything like that for >> the Mat stuff yet. We need to understand the mathematicas better. Pardon my ignorance. What are spinor operations? Why isn't something like a M.*X a matrix operation. It is the Hadamard product. vishy ----- Don't you wish that all the people who sincerely want to help you could agree with each other?

1 0

Fortran 90 interfaces with bfort?
by Barry Smith 19 Aug '05

19 Aug '05

Bill, Do you think it would be possible to enhance bfort to generate Fortran 90 interface definitions? Something like PetscErrorCode MatFactorInfoInitialize(MatFactorInfo info) becomes interface subroutine MatFactorInfoInitialize(info,ierr) MatFactorInfo :: info(11) PetscErrorCode :: ierr end subroutine end interface I think this would help with Fortran programmers writing correct codes much quicker. Thanks, Barry

1 0