Can I use MatSetBlockSize() for MPIAIJ

[Ryan Yan]

Hi All,
I have a question as follows:

In order to use MatSetValuesBlocked() for a MPIAIJ matrix. I need to call
MatSetBlockSize() when I create the matrix.

so I did the following. Here the blocksize = 5;


Mat *A;
....
MatCreate(MPI_COMM_WORLD,A);
MatSetSizes(*A,m*blocksize,n*blocksize,M*blocksize,N*blocksize);
MatSetType(*A,MATMPIAIJ);
MatSetBlockSize(*A,blocksize);
ierr=MatMPIAIJSetPreallocation(*A,0,ourlens_ptws,0,offlens_ptws);
CHKERRQ(ierr);

ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD,"the bs BEFORE is %d\n", bs);
MatGetBlockSize(*A,&bs);

PetscPrintf(PETSC_COMM_WORLD,"the bs is %d\n", bs);
PetscPrintf(PETSC_COMM_WORLD,"the blocksize is %d\n", blocksize);
...

The output I get is:

the bs BEFORE is 0
the bs is 1
the blocksize is 5


It seems like the Mat A does not absorb the information blocksize=5 at all.
How should I make the function-call sequence correct, if I want to set a
blocksize for the MPIAIJ.

Thanks for any suggestions in advance,

Yan

[Ryan Yan]

Sorry a typo:
MatCreate(MPI_COMM_WORLD,*A);

Post by Ryan Yan
Hi All,
In order to use MatSetValuesBlocked() for a MPIAIJ matrix. I need to call
MatSetBlockSize() when I create the matrix.
so I did the following. Here the blocksize = 5;
Mat *A;
....
MatCreate(MPI_COMM_WORLD,A);
MatSetSizes(*A,m*blocksize,n*blocksize,M*blocksize,N*blocksize);
MatSetType(*A,MATMPIAIJ);
MatSetBlockSize(*A,blocksize);
ierr=MatMPIAIJSetPreallocation(*A,0,ourlens_ptws,0,offlens_ptws);
CHKERRQ(ierr);
ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD,"the bs BEFORE is %d\n", bs);
MatGetBlockSize(*A,&bs);
PetscPrintf(PETSC_COMM_WORLD,"the bs is %d\n", bs);
PetscPrintf(PETSC_COMM_WORLD,"the blocksize is %d\n", blocksize);
...
the bs BEFORE is 0
the bs is 1
the blocksize is 5
It seems like the Mat A does not absorb the information blocksize=5 at all.
How should I make the function-call sequence correct, if I want to set a
blocksize for the MPIAIJ.
Thanks for any suggestions in advance,
Yan

[j***@ascomp.ch]

Hello,

Is there an equivalent way to allocating array pointer for creating
vectors (or vector pointers)?


Jarunan

[Matthew Knepley]

Is this what you want?

http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Vec/VecCreateMPIWithArray.html

Matt

Post by j***@ascomp.ch
Hello,
Is there an equivalent way to allocating array pointer for creating vectors
(or vector pointers)?
Jarunan

[Jed Brown]

Post by j***@ascomp.ch
Hello,
Is there an equivalent way to allocating array pointer for creating
vectors (or vector pointers)?

What do you want to do? Maybe you're looking for one of the
VecCreateXXWithArray() variants.

Jed

[j***@ascomp.ch]

I am solving multi-level grid (similar to multi grid but not the
same). In each iteration, each level is solved separately but
solutions are mapped to eachother. Each level has different size of
matrix and vector. And each test case has different numbers of grid
level.

I have a difficulty to create vectors and matrices for each level,
preparing for the computation, as I do not want to create and destroy
them in every iteration. I am thinking of something similar to array
pointer (the code is in fortran) e.g.,

Type(real_array), Dimension(:), allocatable:: pointername
Allocate(pointername(level_numbers))

do i = 1, level_numbers
allocate(pointername(i)%p(size))
enddo


Is it possible to create pointer to vectors?


Thank you
Jarunan

Post by j***@ascomp.ch
Hello,
Is there an equivalent way to allocating array pointer for creating
vectors (or vector pointers)?
Jarunan

[Matthew Knepley]

Yes, Vec is just a regular type.

Matt

I am solving multi-level grid (similar to multi grid but not the same). In
each iteration, each level is solved separately but solutions are mapped to
eachother. Each level has different size of matrix and vector. And each test
case has different numbers of grid level.
I have a difficulty to create vectors and matrices for each level,
preparing for the computation, as I do not want to create and destroy them
in every iteration. I am thinking of something similar to array pointer (the
code is in fortran) e.g.,
Type(real_array), Dimension(:), allocatable:: pointername
Allocate(pointername(level_numbers))
do i = 1, level_numbers
allocate(pointername(i)%p(size))
enddo
Is it possible to create pointer to vectors?
Thank you
Jarunan

Post by j***@ascomp.ch
Hello,
Is there an equivalent way to allocating array pointer for creating
vectors (or vector pointers)?
Jarunan

--
Jarunan Panyasantisuk
Development Engineer
ASCOMP GmbH, Technoparkstr. 1
CH-8005 Zurich, Switzerland
Phone : [phone]
Fax : [phone]
www.ascomp.ch

[Barry Smith]

Yes, for example

integer localn(level_numbers)
Vec myvecs(level_numbers)

for i=1,level_numbers

VecCreate(PETSC_COMM_WORLD,localn(i),PETSC_DETERMINE,myvecs(i),ierr)
endo

Of course, myvecs() can also may made allocatable and you can set
at run time the number of levels.

Barry

Post by j***@ascomp.ch
I am solving multi-level grid (similar to multi grid but not the
same). In each iteration, each level is solved separately but
solutions are mapped to eachother. Each level has different size of
matrix and vector. And each test case has different numbers of grid
level.
I have a difficulty to create vectors and matrices for each level,
preparing for the computation, as I do not want to create and
destroy them in every iteration. I am thinking of something similar
to array pointer (the code is in fortran) e.g.,
Type(real_array), Dimension(:), allocatable:: pointername
Allocate(pointername(level_numbers))
do i = 1, level_numbers
allocate(pointername(i)%p(size))
enddo
Is it possible to create pointer to vectors?
Thank you
Jarunan

Post by j***@ascomp.ch
Hello,
Is there an equivalent way to allocating array pointer for creating
vectors (or vector pointers)?
Jarunan

--
Jarunan Panyasantisuk
Development Engineer
ASCOMP GmbH, Technoparkstr. 1
CH-8005 Zurich, Switzerland
Phone : [phone]
Fax : [phone]
www.ascomp.ch

[Jed Brown]

Post by Ryan Yan
Hi All,
In order to use MatSetValuesBlocked() for a MPIAIJ matrix. I need to
call MatSetBlockSize() when I create the matrix.

Call MatSetBlockSize *after* preallocation.

Jed

[Ryan Yan]

Hi Jed,
Thanks,
So the following is very confusing:

http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatSetValuesBlocked.html

Notes The m and n count the NUMBER of blocks in the row direction and column
direction, NOT the total number of rows/columns; for example, if the block
size<http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Sys/size.html#size>is
2 and you are passing in values for rows 2,3,4,5 then m would be 2
(not
4). The values in idxm would be 1 2; that is the first index for each block
divided by the block
size<http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Sys/size.html#size>.


Note that you must call
MatSetBlockSize<http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatSetBlockSize.html#MatSetBlockSize>()
when constructing this matrix (and before preallocating it).........

Post by Jed Brown

Post by Ryan Yan
Hi All,
In order to use MatSetValuesBlocked() for a MPIAIJ matrix. I need to
call MatSetBlockSize() when I create the matrix.

Call MatSetBlockSize *after* preallocation.
Jed

[Ryan Yan]

It works.
Thanks agian.

Yan

Post by Ryan Yan
Hi Jed,
Thanks,
http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatSetValuesBlocked.html
Notes The m and n count the NUMBER of blocks in the row direction and
column direction, NOT the total number of rows/columns; for example, if the
block size<http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Sys/size.html#size>is 2 and you are passing in values for rows 2,3,4,5 then m would be 2 (not
4). The values in idxm would be 1 2; that is the first index for each block
divided by the block size<http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Sys/size.html#size>.
Note that you must call MatSetBlockSize<http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatSetBlockSize.html#MatSetBlockSize>()
when constructing this matrix (and before preallocating it).........

Post by Jed Brown

Post by Ryan Yan
Hi All,
In order to use MatSetValuesBlocked() for a MPIAIJ matrix. I need to
call MatSetBlockSize() when I create the matrix.

Call MatSetBlockSize *after* preallocation.
Jed

[Jed Brown]

Post by Ryan Yan
Note that you must call MatSetBlockSize
<http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatSetBlockSize.html#MatSetBlockSize>()
when constructing this matrix (and before preallocating it).........

Indeed, thanks for pointing it out. I have fixed the documentation in
petsc-dev and also made MatSetBlockSize() work for BAIJ (it just checks
that the block size agrees with the way the matrix was allocated).

Jed