Re: [petsc-users] Question about MATMPIAIJ and column indexes ordering

2024-01-08 Thread Barry Smith

   Added clarification to the man pages in 
https://gitlab.com/petsc/petsc/-/merge_requests/7170


> On Jan 8, 2024, at 4:31 AM, Deuse, Mathieu via petsc-users 
>  wrote:
> 
> Hello,
>  
> I have a piece of code which generates a matrix in CSR format, but the 
> without sorting the column indexes in increasing order within each row. This 
> seems not to be 100% compatible with the MATMPIAIJ format: the documentation 
> of MatCreateMPIAIJWithArrays indeed mentions 'row-major ordering'.
>  
> For example, consider the 2x2 matrix (1 2; 3 4), which in my code could be 
> stored as i=[0, 2, 4], j=[1, 0, 0, 1], v=[2, 1, 3, 4]. I can generate the 
> matrix as follows (on 1 proc): MatCreateMPIAIJWithArrays(PETSC_COMM_SELF, 2, 
> 2, 2, 2, i, j, v, ). This appears to work fine, and I can then use the 
> matrix in a KSP for example. However, if I try to update the entry values 
> (same order and values v=[2, 1, 3, 4]) with MatUpdateMPIAIJWithArray(matrix, 
> v), it seems that PETSc does not memorize the order of the column indexes and 
> the matrix that I get now is (2 1; 3 4). I get the same result with 
> MatUpdateMPIAIJWithArrays(matrix, 2, 2, 2, 2, i, j, v). On the other hand, if 
> the column indexes are sorted within each row (i=[0, 2, 4], j=[0, 1, 0, 1], 
> v=[1, 2, 3, 4]), then it works fine. I have attached a minimal working 
> example (C++).
>  
> Can I safely rely on MatCreateMPIAIJWithArrays working fine with unsorted 
> column indexes as long as I do not use MatUpdateMPIAIJWithArray(s)? Or should 
> I do the sorting myself before calling MatCreateMPIAIJWithArrays? (or 
> alternatively use another matrix format).
>  
> Thanks in advance for the help.
>  
> Kind regards,
>  
> Mathieu Deuse
> 



Re: [petsc-users] Question about MATMPIAIJ and column indexes ordering

2024-01-08 Thread Barry Smith


> On Jan 8, 2024, at 4:31 AM, Deuse, Mathieu via petsc-users 
>  wrote:
> 
> Hello,
>  
> I have a piece of code which generates a matrix in CSR format, but the 
> without sorting the column indexes in increasing order within each row. This 
> seems not to be 100% compatible with the MATMPIAIJ format: the documentation 
> of MatCreateMPIAIJWithArrays indeed mentions 'row-major ordering'.
>  
> For example, consider the 2x2 matrix (1 2; 3 4), which in my code could be 
> stored as i=[0, 2, 4], j=[1, 0, 0, 1], v=[2, 1, 3, 4]. I can generate the 
> matrix as follows (on 1 proc): MatCreateMPIAIJWithArrays(PETSC_COMM_SELF, 2, 
> 2, 2, 2, i, j, v, ). This appears to work fine, and I can then use the 
> matrix in a KSP for example. However, if I try to update the entry values 
> (same order and values v=[2, 1, 3, 4]) with MatUpdateMPIAIJWithArray(matrix, 
> v), it seems that PETSc does not memorize the order of the column indexes and 
> the matrix that I get now is (2 1; 3 4). I get the same result with 
> MatUpdateMPIAIJWithArrays(matrix, 2, 2, 2, 2, i, j, v). On the other hand, if 
> the column indexes are sorted within each row (i=[0, 2, 4], j=[0, 1, 0, 1], 
> v=[1, 2, 3, 4]), then it works fine. I have attached a minimal working 
> example (C++).
>  
> Can I safely rely on MatCreateMPIAIJWithArrays working fine with unsorted 
> column indexes as long as I do not use MatUpdateMPIAIJWithArray(s)?

   Yes, this is correct. The column indices do not need to be sorted if you 
never call MatUpdateMPIAIJWithArray().



> Or should I do the sorting myself before calling MatCreateMPIAIJWithArrays? 
> (or alternatively use another matrix format).
>  
> Thanks in advance for the help.
>  
> Kind regards,
>  
> Mathieu Deuse
> 



[petsc-users] Question about MATMPIAIJ and column indexes ordering

2024-01-08 Thread Deuse, Mathieu via petsc-users
Hello,

I have a piece of code which generates a matrix in CSR format, but the without 
sorting the column indexes in increasing order within each row. This seems not 
to be 100% compatible with the MATMPIAIJ format: the documentation of 
MatCreateMPIAIJWithArrays indeed mentions 'row-major ordering'.

For example, consider the 2x2 matrix (1 2; 3 4), which in my code could be 
stored as i=[0, 2, 4], j=[1, 0, 0, 1], v=[2, 1, 3, 4]. I can generate the 
matrix as follows (on 1 proc): MatCreateMPIAIJWithArrays(PETSC_COMM_SELF, 2, 2, 
2, 2, i, j, v, ). This appears to work fine, and I can then use the 
matrix in a KSP for example. However, if I try to update the entry values (same 
order and values v=[2, 1, 3, 4]) with MatUpdateMPIAIJWithArray(matrix, v), it 
seems that PETSc does not memorize the order of the column indexes and the 
matrix that I get now is (2 1; 3 4). I get the same result with 
MatUpdateMPIAIJWithArrays(matrix, 2, 2, 2, 2, i, j, v). On the other hand, if 
the column indexes are sorted within each row (i=[0, 2, 4], j=[0, 1, 0, 1], 
v=[1, 2, 3, 4]), then it works fine. I have attached a minimal working example 
(C++).

Can I safely rely on MatCreateMPIAIJWithArrays working fine with unsorted 
column indexes as long as I do not use MatUpdateMPIAIJWithArray(s)? Or should I 
do the sorting myself before calling MatCreateMPIAIJWithArrays? (or 
alternatively use another matrix format).

Thanks in advance for the help.

Kind regards,

Mathieu Deuse
#include 

#include "petsc.h"

int main(int argc, char** argv)
{
PetscCall(PetscInitialize(, , nullptr, nullptr));

PetscIntm   = 2;
PetscIntn   = 2;
PetscInti[] = {0, 2, 4};
PetscIntj[] = {1, 0, 0, 1};
PetscScalar v[] = {2, 1, 3, 4};

Mat mat;
PetscCall(MatCreateMPIAIJWithArrays(PETSC_COMM_SELF, m, n, m, n, i, j, v, 
));

std::cout << "=== Before update ===" << std::endl;
for (PetscInt row = 0; row < m; ++row)
{
for (PetscInt col = 0; col < n; ++col)
{
PetscScalar x;
PetscCall(MatGetValue(mat, row, col, ));
std::cout << ' ' << x;
}
std::cout << std::endl;
}

PetscCall(MatUpdateMPIAIJWithArray(mat, v));

std::cout << "=== After update ===" << std::endl;
for (PetscInt row = 0; row < m; ++row)
{
for (PetscInt col = 0; col < n; ++col)
{
PetscScalar x;
PetscCall(MatGetValue(mat, row, col, ));
std::cout << ' ' << x;
}
std::cout << std::endl;
}

PetscCall(MatDestroy());

PetscCall(PetscFinalize());
}