Re: [petsc-users] Seg fault in gdb but program runs

2022-12-12 Thread Yuyun Yang
Ok I’ll check that, thanks for taking a look! By the way, when I reduce the 
domain size this error doesn’t appear anymore, so I don’t know whether gdb just 
cannot handle the memory, and start to cut things off which is causing the seg 
fault.

From: Matthew Knepley 
Date: Tuesday, December 13, 2022 at 2:49 PM
To: Yuyun Yang 
Cc: petsc-users 
Subject: Re: [petsc-users] Seg fault in gdb but program runs
On Tue, Dec 13, 2022 at 1:14 AM Yuyun Yang 
mailto:yyan...@alumni.stanford.edu>> wrote:
Here is the error message:

Program received signal SIGSEGV, Segmentation fault.
0x555e73b7 in kronConvert (left=..., right=...,
mat=@0x55927e10: 0x57791bb0, diag=5, offDiag=0)
at /home/yuyun/scycle-2/source/spmat.cpp:265
265  kronConvert_symbolic(left,right,mat,d_nnz,o_nnz);

d_nnz and o_nnz are pointers, and they are supposed to hold arrays of the 
number of nonzero in each row,
You seem to be passing integers.

  Thanks,

 Matt

On Tue, Dec 13, 2022 at 12:41 PM Matthew Knepley 
mailto:knep...@gmail.com>> wrote:
On Mon, Dec 12, 2022 at 9:56 PM Yuyun Yang 
mailto:yyan...@alumni.stanford.edu>> wrote:
Hello team,

I’m debugging my code using gdb. The program runs just fine if I don’t debug 
it, but when I use gdb, it seg faults at a place where it never experienced any 
seg fault when I debugged it 1-2 years ago. I wonder if this might be caused by 
the PETSc version change?

The only PETSc calls are the MatGetOwnershipRange() calls, which have not 
changed, so I think this is unlikely.

Or something wrong with gdb itself? I’ve included the code block that is 
problematic for you to take a look at what might be wrong – seg fault happens 
when this function is called. For context, Spmat is a class of sparse matrices 
in the code:

What is the debugger output?

  Thanks,

Matt


// calculate the exact nonzero structure which results from the kronecker outer 
product of left and right

// d_nnz = diagonal nonzero structure, o_nnz = off-diagonal nonzero structure

void kronConvert_symbolic(const Spmat , const Spmat , Mat , 
PetscInt* d_nnz, PetscInt* o_nnz)

{

  size_t rightRowSize = right.size(1);

  size_t rightColSize = right.size(2);



  PetscInt Istart,Iend; // rows owned by current processor

  PetscInt Jstart,Jend; // cols owned by current processor



  // allocate space for mat

  MatGetOwnershipRange(mat,,);

  MatGetOwnershipRangeColumn(mat,,);

  PetscInt m = Iend - Istart;



  for (int ii=0; iisecond).begin(); JjL!=(IiL->second).end(); JjL++) {
  rowL = IiL->first;
  colL = JjL->first;
  valL = JjL->second;
  if (valL==0) { continue; }

  // loop over all values in right
  for (IiR=right._mat.begin(); IiR!=right._mat.end(); IiR++) {
for (JjR=(IiR->second).begin(); JjR!=(IiR->second).end(); JjR++) {
  rowR = IiR->first;
  colR = JjR->first;
  valR = JjR->second;

  // the new values and coordinates for the product matrix
  val = valL*valR;
  row = rowL*rightRowSize + rowR;
  col = colL*rightColSize + colR;

  PetscInt ii = row - Istart; // array index for d_nnz and o_nnz
  if (val!=0 && row >= Istart && row < Iend && col >= Jstart && col < 
Jend) { d_nnz[ii]++; \
}
  if ( (val!=0 && row >= Istart && row < Iend) && (col < Jstart || col 
>= Jend) ) { o_nnz[i\
i]++; }
}
  }
}
  }
}



Thank you,
Yuyun


--
What most experimenters take for granted before they begin their experiments is 
infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/


--
What most experimenters take for granted before they begin their experiments is 
infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/


Re: [petsc-users] Seg fault in gdb but program runs

2022-12-12 Thread Matthew Knepley
On Tue, Dec 13, 2022 at 1:14 AM Yuyun Yang 
wrote:

> Here is the error message:
>
> Program received signal SIGSEGV, Segmentation fault.
> 0x555e73b7 in kronConvert (left=..., right=...,
> mat=@0x55927e10: 0x57791bb0, diag=5, offDiag=0)
> at /home/yuyun/scycle-2/source/spmat.cpp:265
> 265  kronConvert_symbolic(left,right,mat,d_nnz,o_nnz);
>

d_nnz and o_nnz are pointers, and they are supposed to hold arrays of the
number of nonzero in each row,
You seem to be passing integers.

  Thanks,

 Matt


> On Tue, Dec 13, 2022 at 12:41 PM Matthew Knepley 
> wrote:
>
>> On Mon, Dec 12, 2022 at 9:56 PM Yuyun Yang 
>> wrote:
>>
>>> Hello team,
>>>
>>>
>>>
>>> I’m debugging my code using gdb. The program runs just fine if I don’t
>>> debug it, but when I use gdb, it seg faults at a place where it never
>>> experienced any seg fault when I debugged it 1-2 years ago. I wonder if
>>> this might be caused by the PETSc version change?
>>>
>>
>> The only PETSc calls are the MatGetOwnershipRange() calls, which have not
>> changed, so I think this is unlikely.
>>
>>
>>> Or something wrong with gdb itself? I’ve included the code block that is
>>> problematic for you to take a look at what might be wrong – seg fault
>>> happens when this function is called. For context, Spmat is a class of
>>> sparse matrices in the code:
>>>
>>
>> What is the debugger output?
>>
>>   Thanks,
>>
>> Matt
>>
>>
>>> // calculate the exact nonzero structure which results from the
>>> kronecker outer product of left and right
>>>
>>>
>>> // d_nnz = diagonal nonzero structure, o_nnz = off-diagonal nonzero
>>> structure
>>>
>>> void kronConvert_symbolic(const Spmat , const Spmat , Mat &
>>> mat, PetscInt* d_nnz, PetscInt* o_nnz)
>>>
>>>
>>> {
>>>
>>>
>>>   size_t rightRowSize = right.size(1);
>>>
>>>
>>>   size_t rightColSize = right.size(2);
>>>
>>>
>>>
>>>
>>>
>>>   PetscInt Istart,Iend; // rows owned by current processor
>>>
>>>
>>>   PetscInt Jstart,Jend; // cols owned by current processor
>>>
>>>
>>>
>>>
>>>
>>>   // allocate space for mat
>>>
>>>
>>>   MatGetOwnershipRange(mat,,);
>>>
>>>
>>>   MatGetOwnershipRangeColumn(mat,,);
>>>
>>>
>>>   PetscInt m = Iend - Istart;
>>>
>>>
>>>
>>>
>>>
>>>   for (int ii=0; ii>>
>>>
>>>   for (int ii=0; ii>>
>>>
>>>
>>>
>>>
>>>   // iterate over only nnz entries
>>>
>>>
>>>   Spmat::const_row_iter IiL,IiR;
>>>
>>>
>>>   Spmat::const_col_iter JjL,JjR;
>>>
>>>
>>>   double valL=0, valR=0, val=0;
>>>
>>>
>>>   PetscInt row,col;
>>>
>>>
>>>   size_t rowL,colL,rowR,colR;
>>>
>>>
>>>
>>>
>>>   // loop over all values in left
>>>
>>>
>>>   for (IiL=left._mat.begin(); IiL!=left._mat.end(); IiL++) {
>>>
>>>
>>> for (JjL=(IiL->second).begin(); JjL!=(IiL->second).end(); JjL++) {
>>>
>>>
>>>   rowL = IiL->first;
>>>
>>>
>>>   colL = JjL->first;
>>>
>>>
>>>   valL = JjL->second;
>>>
>>>
>>>   if (valL==0) { continue; }
>>>
>>>
>>>
>>>
>>>
>>>   // loop over all values in right
>>>
>>>
>>>   for (IiR=right._mat.begin(); IiR!=right._mat.end(); IiR++) {
>>>
>>>
>>> for (JjR=(IiR->second).begin(); JjR!=(IiR->second).end();
>>> JjR++) {
>>>
>>>   rowR = IiR->first;
>>>
>>>
>>>   colR = JjR->first;
>>>
>>>
>>>   valR = JjR->second;
>>>
>>>
>>>
>>>
>>>
>>>   // the new values and coordinates for the product matrix
>>>
>>>
>>>   val = valL*valR;
>>>
>>>
>>>   row = rowL*rightRowSize + rowR;
>>>
>>>
>>>   col = colL*rightColSize + colR;
>>>
>>>
>>>
>>>
>>>
>>>   PetscInt ii = row - Istart; // array index for d_nnz and
>>> o_nnz
>>>
>>>   if (val!=0 && row >= Istart && row < Iend && col >= Jstart &&
>>> col < Jend) { d_nnz[ii]++; \
>>>
>>> }
>>>
>>>
>>>   if ( (val!=0 && row >= Istart && row < Iend) && (col < Jstart
>>> || col >= Jend) ) { o_nnz[i\
>>>
>>> i]++; }
>>>
>>>
>>> }
>>>
>>>
>>>   }
>>>
>>>
>>> }
>>>
>>>
>>>   }
>>>
>>>
>>> }
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> Thank you,
>>>
>>> Yuyun
>>>
>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> 
>>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 


Re: [petsc-users] Seg fault in gdb but program runs

2022-12-12 Thread Yuyun Yang
Here is the error message:

Program received signal SIGSEGV, Segmentation fault.
0x555e73b7 in kronConvert (left=..., right=...,
mat=@0x55927e10: 0x57791bb0, diag=5, offDiag=0)
at /home/yuyun/scycle-2/source/spmat.cpp:265
265  kronConvert_symbolic(left,right,mat,d_nnz,o_nnz);

On Tue, Dec 13, 2022 at 12:41 PM Matthew Knepley  wrote:

> On Mon, Dec 12, 2022 at 9:56 PM Yuyun Yang 
> wrote:
>
>> Hello team,
>>
>>
>>
>> I’m debugging my code using gdb. The program runs just fine if I don’t
>> debug it, but when I use gdb, it seg faults at a place where it never
>> experienced any seg fault when I debugged it 1-2 years ago. I wonder if
>> this might be caused by the PETSc version change?
>>
>
> The only PETSc calls are the MatGetOwnershipRange() calls, which have not
> changed, so I think this is unlikely.
>
>
>> Or something wrong with gdb itself? I’ve included the code block that is
>> problematic for you to take a look at what might be wrong – seg fault
>> happens when this function is called. For context, Spmat is a class of
>> sparse matrices in the code:
>>
>
> What is the debugger output?
>
>   Thanks,
>
> Matt
>
>
>> // calculate the exact nonzero structure which results from the kronecker
>> outer product of left and right
>>
>>
>> // d_nnz = diagonal nonzero structure, o_nnz = off-diagonal nonzero
>> structure
>>
>> void kronConvert_symbolic(const Spmat , const Spmat , Mat ,
>> PetscInt* d_nnz, PetscInt* o_nnz)
>>
>>
>> {
>>
>>
>>   size_t rightRowSize = right.size(1);
>>
>>
>>   size_t rightColSize = right.size(2);
>>
>>
>>
>>
>>
>>   PetscInt Istart,Iend; // rows owned by current processor
>>
>>
>>   PetscInt Jstart,Jend; // cols owned by current processor
>>
>>
>>
>>
>>
>>   // allocate space for mat
>>
>>
>>   MatGetOwnershipRange(mat,,);
>>
>>
>>   MatGetOwnershipRangeColumn(mat,,);
>>
>>
>>   PetscInt m = Iend - Istart;
>>
>>
>>
>>
>>
>>   for (int ii=0; ii>
>>
>>   for (int ii=0; ii>
>>
>>
>>
>>
>>   // iterate over only nnz entries
>>
>>
>>   Spmat::const_row_iter IiL,IiR;
>>
>>
>>   Spmat::const_col_iter JjL,JjR;
>>
>>
>>   double valL=0, valR=0, val=0;
>>
>>
>>   PetscInt row,col;
>>
>>
>>   size_t rowL,colL,rowR,colR;
>>
>>
>>
>>
>>   // loop over all values in left
>>
>>
>>   for (IiL=left._mat.begin(); IiL!=left._mat.end(); IiL++) {
>>
>>
>> for (JjL=(IiL->second).begin(); JjL!=(IiL->second).end(); JjL++) {
>>
>>
>>   rowL = IiL->first;
>>
>>
>>   colL = JjL->first;
>>
>>
>>   valL = JjL->second;
>>
>>
>>   if (valL==0) { continue; }
>>
>>
>>
>>
>>
>>   // loop over all values in right
>>
>>
>>   for (IiR=right._mat.begin(); IiR!=right._mat.end(); IiR++) {
>>
>>
>> for (JjR=(IiR->second).begin(); JjR!=(IiR->second).end(); JjR++)
>> {
>>
>>   rowR = IiR->first;
>>
>>
>>   colR = JjR->first;
>>
>>
>>   valR = JjR->second;
>>
>>
>>
>>
>>
>>   // the new values and coordinates for the product matrix
>>
>>
>>   val = valL*valR;
>>
>>
>>   row = rowL*rightRowSize + rowR;
>>
>>
>>   col = colL*rightColSize + colR;
>>
>>
>>
>>
>>
>>   PetscInt ii = row - Istart; // array index for d_nnz and o_nnz
>>
>>
>>   if (val!=0 && row >= Istart && row < Iend && col >= Jstart &&
>> col < Jend) { d_nnz[ii]++; \
>>
>> }
>>
>>
>>   if ( (val!=0 && row >= Istart && row < Iend) && (col < Jstart
>> || col >= Jend) ) { o_nnz[i\
>>
>> i]++; }
>>
>>
>> }
>>
>>
>>   }
>>
>>
>> }
>>
>>
>>   }
>>
>>
>> }
>>
>>
>>
>>
>>
>>
>>
>> Thank you,
>>
>> Yuyun
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> 
>


Re: [petsc-users] parallelize matrix assembly process

2022-12-12 Thread Matthew Knepley
On Mon, Dec 12, 2022 at 11:54 PM 김성익  wrote:

>
> With the following example
> https://gitlab.com/petsc/petsc/-/blob/main/src/mat/tests/ex230.c
>
> There are some questions about MatPreallocator.
>
> 1. In parallel run, all the MPI ranks should do the same preallocator
> procedure?
>

In parallel, each process adds its own entries, just as you would in the
matrix assembly.


> 2. In ex230.c, the difference between ex1 of ex230.c and ex2 of ex230.c is
> the block.
> Developers want  to show using block is more efficient method than just
> using matsetvalues?
>

It can be.

  Thanks,

Matt


> Thanks,
> Hyung Kim
>
> 2022년 12월 13일 (화) 오전 1:43, Junchao Zhang 님이 작성:
>
>> Since you run with multiple ranks, you should use matrix type mpiaij and
>> MatMPIAIJSetPreallocation. If preallocation is difficult to estimate, you
>> can use MatPreallocator, see an example at
>> https://gitlab.com/petsc/petsc/-/blob/main/src/mat/tests/ex230.c
>>
>> --Junchao Zhang
>>
>>
>> On Mon, Dec 12, 2022 at 5:16 AM 김성익  wrote:
>>
>>> Hello,
>>>
>>>
>>> I need some keyword or some examples for parallelizing matrix assemble
>>> process.
>>>
>>> My current state is as below.
>>> - Finite element analysis code for Structural mechanics.
>>> - problem size : 3D solid hexa element (number of elements : 125,000),
>>> number of degree of freedom : 397,953
>>> - Matrix type : seqaij, matrix set preallocation by using
>>> MatSeqAIJSetPreallocation
>>> - Matrix assemble time by using 1 core : 120 sec
>>>for (int i=0; i<125000; i++) {
>>> ~~ element matrix calculation}
>>>matassemblybegin
>>>matassemblyend
>>> - Matrix assemble time by using 8 core : 70,234sec
>>>   int start, end;
>>>   VecGetOwnershipRange( element_vec, , );
>>>   for (int i=start; i>>~~ element matrix calculation
>>>matassemblybegin
>>>matassemblyend
>>>
>>>
>>> As you see the state, the parallel case spent a lot of time than
>>> sequential case..
>>> How can I speed up in this case?
>>> Can I get some keyword or examples for parallelizing assembly of matrix
>>> in finite element analysis ?
>>>
>>> Thanks,
>>> Hyung Kim
>>>
>>>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 


Re: [petsc-users] parallelize matrix assembly process

2022-12-12 Thread 김성익
With the following example
https://gitlab.com/petsc/petsc/-/blob/main/src/mat/tests/ex230.c

There are some questions about MatPreallocator.

1. In parallel run, all the MPI ranks should do the same preallocator
procedure?

2. In ex230.c, the difference between ex1 of ex230.c and ex2 of ex230.c is
the block.
Developers want  to show using block is more efficient method than just
using matsetvalues?

Thanks,
Hyung Kim

2022년 12월 13일 (화) 오전 1:43, Junchao Zhang 님이 작성:

> Since you run with multiple ranks, you should use matrix type mpiaij and
> MatMPIAIJSetPreallocation. If preallocation is difficult to estimate, you
> can use MatPreallocator, see an example at
> https://gitlab.com/petsc/petsc/-/blob/main/src/mat/tests/ex230.c
>
> --Junchao Zhang
>
>
> On Mon, Dec 12, 2022 at 5:16 AM 김성익  wrote:
>
>> Hello,
>>
>>
>> I need some keyword or some examples for parallelizing matrix assemble
>> process.
>>
>> My current state is as below.
>> - Finite element analysis code for Structural mechanics.
>> - problem size : 3D solid hexa element (number of elements : 125,000),
>> number of degree of freedom : 397,953
>> - Matrix type : seqaij, matrix set preallocation by using
>> MatSeqAIJSetPreallocation
>> - Matrix assemble time by using 1 core : 120 sec
>>for (int i=0; i<125000; i++) {
>> ~~ element matrix calculation}
>>matassemblybegin
>>matassemblyend
>> - Matrix assemble time by using 8 core : 70,234sec
>>   int start, end;
>>   VecGetOwnershipRange( element_vec, , );
>>   for (int i=start; i>~~ element matrix calculation
>>matassemblybegin
>>matassemblyend
>>
>>
>> As you see the state, the parallel case spent a lot of time than
>> sequential case..
>> How can I speed up in this case?
>> Can I get some keyword or examples for parallelizing assembly of matrix
>> in finite element analysis ?
>>
>> Thanks,
>> Hyung Kim
>>
>>


Re: [petsc-users] realCoords for DOFs

2022-12-12 Thread Matthew Knepley
On Mon, Dec 12, 2022 at 6:06 PM Yann Jobic  wrote:

> Hi all,
>
> I'm trying to get the coords of the dofs of a DMPlex for a PetscFE
> discretization, for orders greater than 1.
>
> I'm struggling to run dm/impls/plex/tutorials/ex8.c
> I've got the following error (with option -view_coord) :
>

You just need to call

  DMGetCoordinatesLocalSetUp()

before the loop. I try to indicate this in the error message(). I did not
call it in the example
because it is only necessary for output. We do this so that output is not
synchronizing.

  Thanks,

Matt


> [0]PETSC ERROR: - Error Message
> --
> [0]PETSC ERROR: Object is in wrong state
> [0]PETSC ERROR: DMGetCoordinatesLocalSetUp() has not been called
> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.18.2, unknown
> [0]PETSC ERROR:
> /home/jobic/projet/fe-utils/petsc/fetools/cmake-build-debug/ex_closure_petsc
>
> on a  named luke by jobic Mon Dec 12 23:34:37 2022
> [0]PETSC ERROR: Configure options
> --prefix=/local/lib/petsc/3.18/p3/gcc/nompi_hdf5 --with-mpi=0
> --with-debugging=1 --with-blacs=1 --download-zlib,--download-p4est
> --download-hdf5=1 --download-triangle=1 --with-single-library=0
> --with-large-file-io=1 --with-shared-libraries=0 -CFLAGS=" -g -O0"
> -CXXFLAGS=" -g -O0" -FFLAGS=" -g -O0" PETSC_ARCH=nompi_gcc_hdf5
> [0]PETSC ERROR: #1 DMGetCoordinatesLocalNoncollective() at
> /home/devel/src_linux/petsc-3.18.0/src/dm/interface/dmcoordinates.c:621
> [0]PETSC ERROR: #2 DMPlexGetCellCoordinates() at
> /home/devel/src_linux/petsc-3.18.0/src/dm/impls/plex/plexgeometry.c:1291
> [0]PETSC ERROR: #3 main() at
> /home/jobic/projet/fe-utils/petsc/fetools/src/ex_closure_petsc.c:86
> [0]PETSC ERROR: PETSc Option Table entries:
> [0]PETSC ERROR: -dm_plex_box_faces 2,2
> [0]PETSC ERROR: -dm_plex_dim 2
> [0]PETSC ERROR: -dm_plex_simplex 0
> [0]PETSC ERROR: -petscspace_degree 1
> [0]PETSC ERROR: -view_coord
> [0]PETSC ERROR: End of Error Message ---send entire
> error message to petsc-ma...@mcs.anl.gov--
>
> Maybe i've done something wrong ?
>
>
> Moreover, i don't quite understand the function DMPlexGetLocalOffsets,
> and how to use it with DMGetCoordinatesLocal. It seems that
> DMGetCoordinatesLocal do not have the coords of the dofs, but only the
> nodes defining the geometry.
>
> I've made some small modifications of ex8.c, but i still have an error :
> [0]PETSC ERROR: - Error Message
> --
> [0]PETSC ERROR: Invalid argument
> [0]PETSC ERROR: Wrong type of object: Parameter # 1
> [0]PETSC ERROR: WARNING! There are option(s) set that were not used!
> Could be the program crashed before they were used or a spelling
> mistake, etc!
> [0]PETSC ERROR: Option left: name:-sol value: vtk:sol.vtu
> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.18.2, unknown
> [0]PETSC ERROR:
> /home/jobic/projet/fe-utils/petsc/fetools/cmake-build-debug/view_coords
> on a  named luke by jobic Mon Dec 12 23:51:05 2022
> [0]PETSC ERROR: Configure options
> --prefix=/local/lib/petsc/3.18/p3/gcc/nompi_hdf5 --with-mpi=0
> --with-debugging=1 --with-blacs=1 --download-zlib,--download-p4est
> --download-hdf5=1 --download-triangle=1 --with-single-library=0
> --with-large-file-io=1 --with-shared-libraries=0 -CFLAGS=" -g -O0"
> -CXXFLAGS=" -g -O0" -FFLAGS=" -g -O0" PETSC_ARCH=nompi_gcc_hdf5
> [0]PETSC ERROR: #1 PetscFEGetHeightSubspace() at
> /home/devel/src_linux/petsc-3.18.0/src/dm/dt/fe/interface/fe.c:1692
> [0]PETSC ERROR: #2 DMPlexGetLocalOffsets() at
> /home/devel/src_linux/petsc-3.18.0/src/dm/impls/plex/plexceed.c:98
> [0]PETSC ERROR: #3 ViewOffsets() at
> /home/jobic/projet/fe-utils/petsc/fetools/src/view_coords.c:28
> [0]PETSC ERROR: #4 main() at
> /home/jobic/projet/fe-utils/petsc/fetools/src/view_coords.c:99
> [0]PETSC ERROR: PETSc Option Table entries:
> [0]PETSC ERROR: -heat_petscspace_degree 2
> [0]PETSC ERROR: -sol vtk:sol.vtu
> [0]PETSC ERROR: End of Error Message ---send entire
> error message to petsc-ma...@mcs.anl.gov--
>
>
> Is dm/impls/plex/tutorials/ex8.c a good example for viewing the coords
> of the dofs of a DMPlex ?
>
>
> Thanks,
>
> Yann
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 


Re: [petsc-users] Seg fault in gdb but program runs

2022-12-12 Thread Matthew Knepley
On Mon, Dec 12, 2022 at 9:56 PM Yuyun Yang 
wrote:

> Hello team,
>
>
>
> I’m debugging my code using gdb. The program runs just fine if I don’t
> debug it, but when I use gdb, it seg faults at a place where it never
> experienced any seg fault when I debugged it 1-2 years ago. I wonder if
> this might be caused by the PETSc version change?
>

The only PETSc calls are the MatGetOwnershipRange() calls, which have not
changed, so I think this is unlikely.


> Or something wrong with gdb itself? I’ve included the code block that is
> problematic for you to take a look at what might be wrong – seg fault
> happens when this function is called. For context, Spmat is a class of
> sparse matrices in the code:
>

What is the debugger output?

  Thanks,

Matt


> // calculate the exact nonzero structure which results from the kronecker
> outer product of left and right
>
>
> // d_nnz = diagonal nonzero structure, o_nnz = off-diagonal nonzero
> structure
>
> void kronConvert_symbolic(const Spmat , const Spmat , Mat ,
> PetscInt* d_nnz, PetscInt* o_nnz)
>
>
> {
>
>
>   size_t rightRowSize = right.size(1);
>
>
>   size_t rightColSize = right.size(2);
>
>
>
>
>
>   PetscInt Istart,Iend; // rows owned by current processor
>
>
>   PetscInt Jstart,Jend; // cols owned by current processor
>
>
>
>
>
>   // allocate space for mat
>
>
>   MatGetOwnershipRange(mat,,);
>
>
>   MatGetOwnershipRangeColumn(mat,,);
>
>
>   PetscInt m = Iend - Istart;
>
>
>
>
>
>   for (int ii=0; ii
>
>   for (int ii=0; ii
>
>
>
>
>   // iterate over only nnz entries
>
>
>   Spmat::const_row_iter IiL,IiR;
>
>
>   Spmat::const_col_iter JjL,JjR;
>
>
>   double valL=0, valR=0, val=0;
>
>
>   PetscInt row,col;
>
>
>   size_t rowL,colL,rowR,colR;
>
>
>
>
>   // loop over all values in left
>
>
>   for (IiL=left._mat.begin(); IiL!=left._mat.end(); IiL++) {
>
>
> for (JjL=(IiL->second).begin(); JjL!=(IiL->second).end(); JjL++) {
>
>
>   rowL = IiL->first;
>
>
>   colL = JjL->first;
>
>
>   valL = JjL->second;
>
>
>   if (valL==0) { continue; }
>
>
>
>
>
>   // loop over all values in right
>
>
>   for (IiR=right._mat.begin(); IiR!=right._mat.end(); IiR++) {
>
>
> for (JjR=(IiR->second).begin(); JjR!=(IiR->second).end(); JjR++)
> {
>
>   rowR = IiR->first;
>
>
>   colR = JjR->first;
>
>
>   valR = JjR->second;
>
>
>
>
>
>   // the new values and coordinates for the product matrix
>
>
>   val = valL*valR;
>
>
>   row = rowL*rightRowSize + rowR;
>
>
>   col = colL*rightColSize + colR;
>
>
>
>
>
>   PetscInt ii = row - Istart; // array index for d_nnz and o_nnz
>
>
>   if (val!=0 && row >= Istart && row < Iend && col >= Jstart &&
> col < Jend) { d_nnz[ii]++; \
>
> }
>
>
>   if ( (val!=0 && row >= Istart && row < Iend) && (col < Jstart
> || col >= Jend) ) { o_nnz[i\
>
> i]++; }
>
>
> }
>
>
>   }
>
>
> }
>
>
>   }
>
>
> }
>
>
>
>
>
>
>
> Thank you,
>
> Yuyun
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 


Re: [petsc-users] parallelize matrix assembly process

2022-12-12 Thread 김성익
Following your comments,
I checked by using '-info'.

As you suspected, most elements being computed on wrong MPI rank.
Also, there are a lot of stashed entries.


Should I divide the domain from the problem define stage?
Or is a proper preallocation sufficient?


[0]  PetscCommDuplicate(): Duplicating a communicator 139687279637472
94370404729840 max tags = 2147483647

[1]  PetscCommDuplicate(): Duplicating a communicator 139620736898016
94891084133376 max tags = 2147483647

[0]  MatSetUp(): Warning not preallocating matrix storage

[1]  PetscCommDuplicate(): Duplicating a communicator 139620736897504
94891083133744 max tags = 2147483647

[0]  PetscCommDuplicate(): Duplicating a communicator 139687279636960
94370403730224 max tags = 2147483647

[1]  PetscCommDuplicate(): Using internal PETSc communicator
139620736897504 94891083133744

[0]  PetscCommDuplicate(): Using internal PETSc communicator
139687279636960 94370403730224

[1]  PetscCommDuplicate(): Using internal PETSc communicator
139620736898016 94891084133376

[0]  PetscCommDuplicate(): Using internal PETSc communicator
139687279637472 94370404729840

[1]  PetscCommDuplicate(): Using internal PETSc communicator
139620736898016 94891084133376

[0]  PetscCommDuplicate(): Using internal PETSc communicator
139687279637472 94370404729840

[1]  PetscCommDuplicate(): Using internal PETSc communicator
139620736898016 94891084133376

[0]  PetscCommDuplicate(): Using internal PETSc communicator
139687279637472 94370404729840

[0]  PetscCommDuplicate(): Using internal PETSc communicator
139687279637472 94370404729840

[1]  PetscCommDuplicate(): Using internal PETSc communicator
139620736898016 94891084133376

[1]  PetscCommDuplicate(): Using internal PETSc communicator
139620736898016 94891084133376

[0]  PetscCommDuplicate(): Using internal PETSc communicator
139687279637472 94370404729840

 TIME0 : 0.00

 TIME0 : 0.00

[0]  VecAssemblyBegin_MPI_BTS(): Stash has 661 entries, uses 8 mallocs.

[0]  VecAssemblyBegin_MPI_BTS(): Block-Stash has 0 entries, uses 0
mallocs.

[0]  VecAssemblyBegin_MPI_BTS(): Stash has 661 entries, uses 5 mallocs.

[0]  VecAssemblyBegin_MPI_BTS(): Block-Stash has 0 entries, uses 0
mallocs.

[0]  MatAssemblyBegin_MPIAIJ(): Stash has 460416 entries, uses 5
mallocs.

[1]  MatAssemblyBegin_MPIAIJ(): Stash has 461184 entries, uses 5
mallocs.

[0]  MatAssemblyEnd_SeqAIJ(): Matrix size: 13892 X 13892; storage
space: 180684 unneeded,987406 used

[0]  MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues()
is 73242

[0]  MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 81

[0]  MatCheckCompressedRow(): Found the ratio (num_zerorows
0)/(num_localrows 13892) < 0.6. Do not use CompressedRow routines.

[0]  MatSeqAIJCheckInode(): Found 4631 nodes of 13892. Limit used: 5.
Using Inode routines

[1]  MatAssemblyEnd_SeqAIJ(): Matrix size: 13891 X 13891; storage
space: 180715 unneeded,987325 used

[1]  MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues()
is 73239

[1]  MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 81

[1]  MatCheckCompressedRow(): Found the ratio (num_zerorows
0)/(num_localrows 13891) < 0.6. Do not use CompressedRow routines.

[1]  MatSeqAIJCheckInode(): Found 4631 nodes of 13891. Limit used: 5.
Using Inode routines

[0]  PetscCommDuplicate(): Using internal PETSc communicator
139687279636960 94370403730224

[0]  PetscCommDuplicate(): Using internal PETSc communicator
139687279636960 94370403730224

[1]  PetscCommDuplicate(): Using internal PETSc communicator
139620736897504 94891083133744

[1]  PetscCommDuplicate(): Using internal PETSc communicator
139620736897504 94891083133744

[0]  MatAssemblyEnd_SeqAIJ(): Matrix size: 13892 X 1390; storage
space: 72491 unneeded,34049 used

[0]  MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues()
is 2472

[0]  MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 40

[0]  MatCheckCompressedRow(): Found the ratio (num_zerorows
12501)/(num_localrows 13892) > 0.6. Use CompressedRow routines.

Assemble Time : 174.079366sec

[0]  PetscCommDuplicate(): Using internal PETSc communicator
139687279636960 94370403730224

[0]  PetscCommDuplicate(): Using internal PETSc communicator
139687279636960 94370403730224

[1]  MatAssemblyEnd_SeqAIJ(): Matrix size: 13891 X 1391; storage
space: 72441 unneeded,34049 used

[1]  MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues()
is 2469

[1]  MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 41

[1]  MatCheckCompressedRow(): Found the ratio (num_zerorows
12501)/(num_localrows 13891) > 0.6. Use CompressedRow routines.

Assemble Time : 174.141234sec

[1]  PetscCommDuplicate(): Using internal PETSc communicator
139620736897504 94891083133744

[1]  PetscCommDuplicate(): Using internal PETSc communicator
139620736897504 94891083133744

[0]  VecAssemblyBegin_MPI_BTS(): Stash has 13891 entries, uses 8
mallocs.

[0]  VecAssemblyBegin_MPI_BTS(): Block-Stash has 0 entries, uses 0
mallocs.

[1]  

[petsc-users] Seg fault in gdb but program runs

2022-12-12 Thread Yuyun Yang
Hello team,

I’m debugging my code using gdb. The program runs just fine if I don’t debug 
it, but when I use gdb, it seg faults at a place where it never experienced any 
seg fault when I debugged it 1-2 years ago. I wonder if this might be caused by 
the PETSc version change? Or something wrong with gdb itself? I’ve included the 
code block that is problematic for you to take a look at what might be wrong – 
seg fault happens when this function is called. For context, Spmat is a class 
of sparse matrices in the code:


// calculate the exact nonzero structure which results from the kronecker outer 
product of left and right

// d_nnz = diagonal nonzero structure, o_nnz = off-diagonal nonzero structure

void kronConvert_symbolic(const Spmat , const Spmat , Mat , 
PetscInt* d_nnz, PetscInt* o_nnz)

{

  size_t rightRowSize = right.size(1);

  size_t rightColSize = right.size(2);



  PetscInt Istart,Iend; // rows owned by current processor

  PetscInt Jstart,Jend; // cols owned by current processor



  // allocate space for mat

  MatGetOwnershipRange(mat,,);

  MatGetOwnershipRangeColumn(mat,,);

  PetscInt m = Iend - Istart;



  for (int ii=0; iisecond).begin(); JjL!=(IiL->second).end(); JjL++) {
  rowL = IiL->first;
  colL = JjL->first;
  valL = JjL->second;
  if (valL==0) { continue; }

  // loop over all values in right
  for (IiR=right._mat.begin(); IiR!=right._mat.end(); IiR++) {
for (JjR=(IiR->second).begin(); JjR!=(IiR->second).end(); JjR++) {
  rowR = IiR->first;
  colR = JjR->first;
  valR = JjR->second;

  // the new values and coordinates for the product matrix
  val = valL*valR;
  row = rowL*rightRowSize + rowR;
  col = colL*rightColSize + colR;

  PetscInt ii = row - Istart; // array index for d_nnz and o_nnz
  if (val!=0 && row >= Istart && row < Iend && col >= Jstart && col < 
Jend) { d_nnz[ii]++; \
}
  if ( (val!=0 && row >= Istart && row < Iend) && (col < Jstart || col 
>= Jend) ) { o_nnz[i\
i]++; }
}
  }
}
  }
}



Thank you,
Yuyun


Re: [petsc-users] realCoords for DOFs

2022-12-12 Thread Mark Adams
PETSc does not store the coordinates for high order elements (ie, the
"midside nodes").
I get them by projecting a f(x) = x, function.
Not sure if that is what you want but I can give you a code snippet if
there are no better ideas.

Mark


On Mon, Dec 12, 2022 at 6:06 PM Yann Jobic  wrote:

> Hi all,
>
> I'm trying to get the coords of the dofs of a DMPlex for a PetscFE
> discretization, for orders greater than 1.
>
> I'm struggling to run dm/impls/plex/tutorials/ex8.c
> I've got the following error (with option -view_coord) :
>
> [0]PETSC ERROR: - Error Message
> --
> [0]PETSC ERROR: Object is in wrong state
> [0]PETSC ERROR: DMGetCoordinatesLocalSetUp() has not been called
> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.18.2, unknown
> [0]PETSC ERROR:
> /home/jobic/projet/fe-utils/petsc/fetools/cmake-build-debug/ex_closure_petsc
>
> on a  named luke by jobic Mon Dec 12 23:34:37 2022
> [0]PETSC ERROR: Configure options
> --prefix=/local/lib/petsc/3.18/p3/gcc/nompi_hdf5 --with-mpi=0
> --with-debugging=1 --with-blacs=1 --download-zlib,--download-p4est
> --download-hdf5=1 --download-triangle=1 --with-single-library=0
> --with-large-file-io=1 --with-shared-libraries=0 -CFLAGS=" -g -O0"
> -CXXFLAGS=" -g -O0" -FFLAGS=" -g -O0" PETSC_ARCH=nompi_gcc_hdf5
> [0]PETSC ERROR: #1 DMGetCoordinatesLocalNoncollective() at
> /home/devel/src_linux/petsc-3.18.0/src/dm/interface/dmcoordinates.c:621
> [0]PETSC ERROR: #2 DMPlexGetCellCoordinates() at
> /home/devel/src_linux/petsc-3.18.0/src/dm/impls/plex/plexgeometry.c:1291
> [0]PETSC ERROR: #3 main() at
> /home/jobic/projet/fe-utils/petsc/fetools/src/ex_closure_petsc.c:86
> [0]PETSC ERROR: PETSc Option Table entries:
> [0]PETSC ERROR: -dm_plex_box_faces 2,2
> [0]PETSC ERROR: -dm_plex_dim 2
> [0]PETSC ERROR: -dm_plex_simplex 0
> [0]PETSC ERROR: -petscspace_degree 1
> [0]PETSC ERROR: -view_coord
> [0]PETSC ERROR: End of Error Message ---send entire
> error message to petsc-ma...@mcs.anl.gov--
>
> Maybe i've done something wrong ?
>
>
> Moreover, i don't quite understand the function DMPlexGetLocalOffsets,
> and how to use it with DMGetCoordinatesLocal. It seems that
> DMGetCoordinatesLocal do not have the coords of the dofs, but only the
> nodes defining the geometry.
>
> I've made some small modifications of ex8.c, but i still have an error :
> [0]PETSC ERROR: - Error Message
> --
> [0]PETSC ERROR: Invalid argument
> [0]PETSC ERROR: Wrong type of object: Parameter # 1
> [0]PETSC ERROR: WARNING! There are option(s) set that were not used!
> Could be the program crashed before they were used or a spelling
> mistake, etc!
> [0]PETSC ERROR: Option left: name:-sol value: vtk:sol.vtu
> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.18.2, unknown
> [0]PETSC ERROR:
> /home/jobic/projet/fe-utils/petsc/fetools/cmake-build-debug/view_coords
> on a  named luke by jobic Mon Dec 12 23:51:05 2022
> [0]PETSC ERROR: Configure options
> --prefix=/local/lib/petsc/3.18/p3/gcc/nompi_hdf5 --with-mpi=0
> --with-debugging=1 --with-blacs=1 --download-zlib,--download-p4est
> --download-hdf5=1 --download-triangle=1 --with-single-library=0
> --with-large-file-io=1 --with-shared-libraries=0 -CFLAGS=" -g -O0"
> -CXXFLAGS=" -g -O0" -FFLAGS=" -g -O0" PETSC_ARCH=nompi_gcc_hdf5
> [0]PETSC ERROR: #1 PetscFEGetHeightSubspace() at
> /home/devel/src_linux/petsc-3.18.0/src/dm/dt/fe/interface/fe.c:1692
> [0]PETSC ERROR: #2 DMPlexGetLocalOffsets() at
> /home/devel/src_linux/petsc-3.18.0/src/dm/impls/plex/plexceed.c:98
> [0]PETSC ERROR: #3 ViewOffsets() at
> /home/jobic/projet/fe-utils/petsc/fetools/src/view_coords.c:28
> [0]PETSC ERROR: #4 main() at
> /home/jobic/projet/fe-utils/petsc/fetools/src/view_coords.c:99
> [0]PETSC ERROR: PETSc Option Table entries:
> [0]PETSC ERROR: -heat_petscspace_degree 2
> [0]PETSC ERROR: -sol vtk:sol.vtu
> [0]PETSC ERROR: End of Error Message ---send entire
> error message to petsc-ma...@mcs.anl.gov--
>
>
> Is dm/impls/plex/tutorials/ex8.c a good example for viewing the coords
> of the dofs of a DMPlex ?
>
>
> Thanks,
>
> Yann
>
>


[petsc-users] realCoords for DOFs

2022-12-12 Thread Yann Jobic

Hi all,

I'm trying to get the coords of the dofs of a DMPlex for a PetscFE 
discretization, for orders greater than 1.


I'm struggling to run dm/impls/plex/tutorials/ex8.c
I've got the following error (with option -view_coord) :

[0]PETSC ERROR: - Error Message 
--

[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: DMGetCoordinatesLocalSetUp() has not been called
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.18.2, unknown
[0]PETSC ERROR: 
/home/jobic/projet/fe-utils/petsc/fetools/cmake-build-debug/ex_closure_petsc 
on a  named luke by jobic Mon Dec 12 23:34:37 2022
[0]PETSC ERROR: Configure options 
--prefix=/local/lib/petsc/3.18/p3/gcc/nompi_hdf5 --with-mpi=0 
--with-debugging=1 --with-blacs=1 --download-zlib,--download-p4est 
--download-hdf5=1 --download-triangle=1 --with-single-library=0 
--with-large-file-io=1 --with-shared-libraries=0 -CFLAGS=" -g -O0" 
-CXXFLAGS=" -g -O0" -FFLAGS=" -g -O0" PETSC_ARCH=nompi_gcc_hdf5
[0]PETSC ERROR: #1 DMGetCoordinatesLocalNoncollective() at 
/home/devel/src_linux/petsc-3.18.0/src/dm/interface/dmcoordinates.c:621
[0]PETSC ERROR: #2 DMPlexGetCellCoordinates() at 
/home/devel/src_linux/petsc-3.18.0/src/dm/impls/plex/plexgeometry.c:1291
[0]PETSC ERROR: #3 main() at 
/home/jobic/projet/fe-utils/petsc/fetools/src/ex_closure_petsc.c:86

[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -dm_plex_box_faces 2,2
[0]PETSC ERROR: -dm_plex_dim 2
[0]PETSC ERROR: -dm_plex_simplex 0
[0]PETSC ERROR: -petscspace_degree 1
[0]PETSC ERROR: -view_coord
[0]PETSC ERROR: End of Error Message ---send entire 
error message to petsc-ma...@mcs.anl.gov--


Maybe i've done something wrong ?


Moreover, i don't quite understand the function DMPlexGetLocalOffsets, 
and how to use it with DMGetCoordinatesLocal. It seems that 
DMGetCoordinatesLocal do not have the coords of the dofs, but only the 
nodes defining the geometry.


I've made some small modifications of ex8.c, but i still have an error :
[0]PETSC ERROR: - Error Message 
--

[0]PETSC ERROR: Invalid argument
[0]PETSC ERROR: Wrong type of object: Parameter # 1
[0]PETSC ERROR: WARNING! There are option(s) set that were not used! 
Could be the program crashed before they were used or a spelling 
mistake, etc!

[0]PETSC ERROR: Option left: name:-sol value: vtk:sol.vtu
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.18.2, unknown
[0]PETSC ERROR: 
/home/jobic/projet/fe-utils/petsc/fetools/cmake-build-debug/view_coords 
on a  named luke by jobic Mon Dec 12 23:51:05 2022
[0]PETSC ERROR: Configure options 
--prefix=/local/lib/petsc/3.18/p3/gcc/nompi_hdf5 --with-mpi=0 
--with-debugging=1 --with-blacs=1 --download-zlib,--download-p4est 
--download-hdf5=1 --download-triangle=1 --with-single-library=0 
--with-large-file-io=1 --with-shared-libraries=0 -CFLAGS=" -g -O0" 
-CXXFLAGS=" -g -O0" -FFLAGS=" -g -O0" PETSC_ARCH=nompi_gcc_hdf5
[0]PETSC ERROR: #1 PetscFEGetHeightSubspace() at 
/home/devel/src_linux/petsc-3.18.0/src/dm/dt/fe/interface/fe.c:1692
[0]PETSC ERROR: #2 DMPlexGetLocalOffsets() at 
/home/devel/src_linux/petsc-3.18.0/src/dm/impls/plex/plexceed.c:98
[0]PETSC ERROR: #3 ViewOffsets() at 
/home/jobic/projet/fe-utils/petsc/fetools/src/view_coords.c:28
[0]PETSC ERROR: #4 main() at 
/home/jobic/projet/fe-utils/petsc/fetools/src/view_coords.c:99

[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -heat_petscspace_degree 2
[0]PETSC ERROR: -sol vtk:sol.vtu
[0]PETSC ERROR: End of Error Message ---send entire 
error message to petsc-ma...@mcs.anl.gov--



Is dm/impls/plex/tutorials/ex8.c a good example for viewing the coords 
of the dofs of a DMPlex ?



Thanks,

Yann

static char help[] = "Test the function \n\n";
/* run with -heat_petscspace_degree 2  -sol vtk:sol.vtu */
/* En fait petsc, pour le heat_petscspace_degree 0 garde le maillage ori, 
traingulaire 3 sommets */
/* Pour les ordres superieurs, il raffine les elements, donc a 2 fois plus 
d'elements par directions */
/* et garde de l'ordre 1 en geometrie   
 */
/* Il faut tester avec HDF5 voir si ce format ne permet pas des elements de 
plus haut degres */

#include "petscdm.h"
#include 
#include 
#include 
#include "petscksp.h"

/* simple gaussian centered at (0.5,0.5) */
static void gaussian(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt 
uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar 
u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], 
const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 
PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar 

Re: [petsc-users] Union of sequential vecs

2022-12-12 Thread Matthew Knepley
On Mon, Dec 12, 2022 at 10:00 AM Karthikeyan Chockalingam - STFC UKRI <
karthikeyan.chockalin...@stfc.ac.uk> wrote:

> Thank you Matt and Blaise.
>
>
>
> I will try out IS (though I have not it before)
>
> (i)  Can IS be of different size, on different
> processors, and still call ISALLGather?
>
Yes.

> (ii)Can IS be passed as row indices to MatZeroRowsColumns?
>

https://petsc.org/main/docs/manualpages/Mat/MatZeroRowsColumnsIS/

  Thanks,

 Matt


> I will look into DMlabels and start a different thread if needed.
>
>
>
> Best,
>
> Karthik.
>
>
>
>
>
>
>
> *From: *Matthew Knepley 
> *Date: *Saturday, 10 December 2022 at 02:20
> *To: *Chockalingam, Karthikeyan (STFC,DL,HC) <
> karthikeyan.chockalin...@stfc.ac.uk>
> *Cc: *Mark Adams , petsc-users@mcs.anl.gov <
> petsc-users@mcs.anl.gov>
> *Subject: *Re: [petsc-users] Union of sequential vecs
>
> On Fri, Dec 9, 2022 at 6:50 PM Karthikeyan Chockalingam - STFC UKRI via
> petsc-users  wrote:
>
> Thank you Mark and Barry.
>
>
>
> @Mark Adams  I follow you for the most part. Shouldn’t R
> be an MPI Vector?
>
>
>
> Here it goes:
>
>
>
> Q[Vec[i]] for all "i" in my local (sequential) "Vec".
>
> Compute number of local nonzeros in Q, call it n.
>
>
>
> //Create a MPI Vector R, with local size n
>
>
>
> VecCreate(PETSC_COMM_WORLD, );
>
> VecSetType(R, VECMPI);
>
> VecSetSizes(R, *n*, PETSC_DECIDE);
>
>
>
> //Populate MPI Vector R with local (sequential) "Vec".
>
>
>
> VecGetOwnershipRange(R, , );
>
> local_size = iend - istart; //The local_size should be ‘*n*’ right?
>
> VecGetArray(R, );
>
>   *for* (i = 0; i < local_size; i++) {
>
> values[i] = Vec[i];
>
>   }
>
> VecRestoreArray(R, );
>
>
>
> //Scatter R to all processors
>
>
>
> VecV_SEQ;
>
> VecScatter ctx;
>
>
>
> VecScatterCreateToAll(R,,_SEQ);
>
>
>
> //Remove duplicates in V_SEQ
>
> How can I use PetscSortRemoveDupsInt to remove duplicates in V_SEQ?
>
>
>
>
>
> Physics behind:
>
> I am reading a parallel mesh, and want to mark all the boundary nodes. I
> use a local (sequential) Vec to store the boundary nodes for each parallel
> partition. Hence, local Vecs can end up with duplicate node index among
> them, which I would like to get rid of when I combine all of them together.
>
>
>
> 1) Blaise is right you should use an IS, not a Vec, to hold node indices.
> His solution is only a few lines, so I would use it.
>
>
>
> 2) I would not recommend doing things this way in the first place. PETSc
> can manage parallel meshes scalably, marking boundaries using DMLabel
> objects.
>
>
>
>   Thanks,
>
>
>
>  Matt
>
>
>
> Best,
>
> Karthik.
>
>
>
>
>
>
>
>
>
> *From: *Mark Adams 
> *Date: *Friday, 9 December 2022 at 21:08
> *To: *Chockalingam, Karthikeyan (STFC,DL,HC) <
> karthikeyan.chockalin...@stfc.ac.uk>
> *Cc: *Barry Smith , petsc-users@mcs.anl.gov <
> petsc-users@mcs.anl.gov>
> *Subject: *Re: [petsc-users] Union of sequential vecs
>
> If your space is pretty compact, eg, (0,12), you could create an MPI
> vector Q of size 13, say, and each processor can add 1.0 to Q[Vec[i]], for
> all "i" in my local "Vec".
>
> Then each processor can count the number of local nonzeros in Q, call it
> n, create a new vector, R, with local size n, then set R[i] = global index
> of the nonzero for each nonzero in Q, i=0:n.
>
> Do some sort of vec-scatter-to-all with R to get what you want.
>
>
>
> Does that work?
>
>
>
> Mark
>
>
>
>
>
> On Fri, Dec 9, 2022 at 3:25 PM Karthikeyan Chockalingam - STFC UKRI via
> petsc-users  wrote:
>
> That is where I am stuck, *I don’t know* who to combine them to get Vec =
> {2,5,7,8,10,11,12}.
>
> I just want them in an MPI vector.
>
>
>
> I finally plan to call VecScatterCreateToAll so that all processor gets a
> copy.
>
>
>
> Thank you.
>
>
>
> Kind regards,
>
> Karthik.
>
>
>
> *From: *Barry Smith 
> *Date: *Friday, 9 December 2022 at 20:04
> *To: *Chockalingam, Karthikeyan (STFC,DL,HC) <
> karthikeyan.chockalin...@stfc.ac.uk>
> *Cc: *petsc-users@mcs.anl.gov 
> *Subject: *Re: [petsc-users] Union of sequential vecs
>
>
>
>   How are you combining them to get Vec = {2,5,7,8,10,11,12}?
>
>
>
>   Do you want the values to remain on the same MPI rank as before, just in
> an MPI vector?
>
>
>
>
>
>
>
> On Dec 9, 2022, at 2:28 PM, Karthikeyan Chockalingam - STFC UKRI via
> petsc-users  wrote:
>
>
>
> Hi,
>
>
>
> I want to take the union of a set of sequential vectors, each living in a
> different processor.
>
>
>
> Say,
>
> Vec_Seq1 = {2,5,7}
>
> Vec_Seq2 = {5,8,10,11}
>
> Vec_Seq3 = {5,2,12}.
>
>
>
> Finally, get the union of all them Vec = {2,5,7,8,10,11,12}.
>
>
>
> I initially wanted to create a parallel vector and insert the (sequential
> vector) values but I do not know, to which index to insert the values to.
> But I do know the total size of Vec (which in this case is 7).
>
>
>
> Any help is much appreciated.
>
>
>
> Kind regards,
>
> Karthik.
>
>
>
>
>
>
>
> This email and any attachments are intended solely for the use of the
> 

Re: [petsc-users] Insert one sparse matrix as a block in another

2022-12-12 Thread Matthew Knepley
On Mon, Dec 12, 2022 at 5:24 PM Jed Brown  wrote:

> The description matches MATNEST (MATCOMPOSITE is for a sum or product of
> matrices) or parallel decompositions. Also consider the assembly style of
> src/snes/tutorials/ex28.c, which can create either a monolithic or block
> (MATNEST) matrix without extra storage or conversion costs.
>

I will just say a few words about ex28. The idea is that if you are already
calling MatSetValues() to assemble your submatrices, then
you can use MatSetValuesLocal() to remap those locations into locations in
the large matrix, using a LocalToGlobalMap. This allows
you to choose either a standard AIJ matrix (which supports factorizations
for example), or a MatNest object that supports fast extraction
of the blocks.

  Thanks,

Matt


> Mark Adams  writes:
>
> > Do you know what kind of solver works well for this problem?
> >
> > You probably want to figure that out first and not worry about
> efficiency.
> >
> > MATCOMPOSITE does what you want but not all solvers will work with it.
> >
> > Where does this problem come from? We have a lot of experience and might
> > know something.
> >
> > Mark
> >
> > On Mon, Dec 12, 2022 at 1:33 PM Peder Jørgensgaard Olesen via
> petsc-users <
> > petsc-users@mcs.anl.gov> wrote:
> >
> >> Hello
> >>
> >>
> >> I have a set of sparse matrices (A1, A2, ...) , and need to generate a
> >> larger matrix B with these as submatrices. I do not know the precise
> sparse
> >> layouts of the A's (only that each row has one or two non-zero values),
> >> and extracting *all* values to copy into B seems incredibly wasteful.
> How
> >> can I make use of the sparsity to solve this efficiently?
> >>
> >>
> >> Thanks,
> >>
> >> Peder
> >>
> >>
> >>
> >> Peder Jørgensgaard Olesen
> >> PhD student
> >> Department of Civil and Mechanical Engineering
> >>
> >> pj...@mek.dtu.dk
> >> Koppels Allé
> >> Building 403, room 105
> >> 2800 Kgs. Lyngby
> >> www.dtu.dk/english
> >>
> >>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 


Re: [petsc-users] Insert one sparse matrix as a block in another

2022-12-12 Thread Jed Brown
The description matches MATNEST (MATCOMPOSITE is for a sum or product of 
matrices) or parallel decompositions. Also consider the assembly style of 
src/snes/tutorials/ex28.c, which can create either a monolithic or block 
(MATNEST) matrix without extra storage or conversion costs.

Mark Adams  writes:

> Do you know what kind of solver works well for this problem?
>
> You probably want to figure that out first and not worry about efficiency.
>
> MATCOMPOSITE does what you want but not all solvers will work with it.
>
> Where does this problem come from? We have a lot of experience and might
> know something.
>
> Mark
>
> On Mon, Dec 12, 2022 at 1:33 PM Peder Jørgensgaard Olesen via petsc-users <
> petsc-users@mcs.anl.gov> wrote:
>
>> Hello
>>
>>
>> I have a set of sparse matrices (A1, A2, ...) , and need to generate a
>> larger matrix B with these as submatrices. I do not know the precise sparse
>> layouts of the A's (only that each row has one or two non-zero values),
>> and extracting *all* values to copy into B seems incredibly wasteful. How
>> can I make use of the sparsity to solve this efficiently?
>>
>>
>> Thanks,
>>
>> Peder
>>
>>
>>
>> Peder Jørgensgaard Olesen
>> PhD student
>> Department of Civil and Mechanical Engineering
>>
>> pj...@mek.dtu.dk
>> Koppels Allé
>> Building 403, room 105
>> 2800 Kgs. Lyngby
>> www.dtu.dk/english
>>
>>


Re: [petsc-users] Insert one sparse matrix as a block in another

2022-12-12 Thread Mark Adams
Do you know what kind of solver works well for this problem?

You probably want to figure that out first and not worry about efficiency.

MATCOMPOSITE does what you want but not all solvers will work with it.

Where does this problem come from? We have a lot of experience and might
know something.

Mark

On Mon, Dec 12, 2022 at 1:33 PM Peder Jørgensgaard Olesen via petsc-users <
petsc-users@mcs.anl.gov> wrote:

> Hello
>
>
> I have a set of sparse matrices (A1, A2, ...) , and need to generate a
> larger matrix B with these as submatrices. I do not know the precise sparse
> layouts of the A's (only that each row has one or two non-zero values),
> and extracting *all* values to copy into B seems incredibly wasteful. How
> can I make use of the sparsity to solve this efficiently?
>
>
> Thanks,
>
> Peder
>
>
>
> Peder Jørgensgaard Olesen
> PhD student
> Department of Civil and Mechanical Engineering
>
> pj...@mek.dtu.dk
> Koppels Allé
> Building 403, room 105
> 2800 Kgs. Lyngby
> www.dtu.dk/english
>
>


[petsc-users] Insert one sparse matrix as a block in another

2022-12-12 Thread Peder Jørgensgaard Olesen via petsc-users
Hello


I have a set of sparse matrices (A1, A2, ...) , and need to generate a larger 
matrix B with these as submatrices. I do not know the precise sparse layouts of 
the A's (only that each row has one or two non-zero values), and extracting all 
values to copy into B seems incredibly wasteful. How can I make use of the 
sparsity to solve this efficiently?


Thanks,

Peder



[http://www.dtu.dk/-/media/DTU_Generelt/Andet/mail-signature-logo.png]

Peder Jørgensgaard Olesen
PhD student
Department of Civil and Mechanical Engineering

pj...@mek.dtu.dk
Koppels Allé
Building 403, room 105
2800 Kgs. Lyngby
www.dtu.dk/english




Re: [petsc-users] parallelize matrix assembly process

2022-12-12 Thread Junchao Zhang
Since you run with multiple ranks, you should use matrix type mpiaij and
MatMPIAIJSetPreallocation. If preallocation is difficult to estimate, you
can use MatPreallocator, see an example at
https://gitlab.com/petsc/petsc/-/blob/main/src/mat/tests/ex230.c

--Junchao Zhang


On Mon, Dec 12, 2022 at 5:16 AM 김성익  wrote:

> Hello,
>
>
> I need some keyword or some examples for parallelizing matrix assemble
> process.
>
> My current state is as below.
> - Finite element analysis code for Structural mechanics.
> - problem size : 3D solid hexa element (number of elements : 125,000),
> number of degree of freedom : 397,953
> - Matrix type : seqaij, matrix set preallocation by using
> MatSeqAIJSetPreallocation
> - Matrix assemble time by using 1 core : 120 sec
>for (int i=0; i<125000; i++) {
> ~~ element matrix calculation}
>matassemblybegin
>matassemblyend
> - Matrix assemble time by using 8 core : 70,234sec
>   int start, end;
>   VecGetOwnershipRange( element_vec, , );
>   for (int i=start; i~~ element matrix calculation
>matassemblybegin
>matassemblyend
>
>
> As you see the state, the parallel case spent a lot of time than
> sequential case..
> How can I speed up in this case?
> Can I get some keyword or examples for parallelizing assembly of matrix in
> finite element analysis ?
>
> Thanks,
> Hyung Kim
>
>


Re: [petsc-users] parallelize matrix assembly process

2022-12-12 Thread Barry Smith


   The problem is possibly due to most elements being computed on "wrong" MPI 
rank and thus requiring almost all the matrix entries to be "stashed" when 
computed and then sent off to the owning MPI rank.  Please send ALL the output 
of a parallel run with -info so we can see how much communication is done in 
the matrix assembly.

  Barry


> On Dec 12, 2022, at 6:16 AM, 김성익  wrote:
> 
> Hello,
> 
> 
> I need some keyword or some examples for parallelizing matrix assemble 
> process.
> 
> My current state is as below.
> - Finite element analysis code for Structural mechanics.
> - problem size : 3D solid hexa element (number of elements : 125,000), number 
> of degree of freedom : 397,953
> - Matrix type : seqaij, matrix set preallocation by using 
> MatSeqAIJSetPreallocation
> - Matrix assemble time by using 1 core : 120 sec
>for (int i=0; i<125000; i++) {
> ~~ element matrix calculation}
>matassemblybegin
>matassemblyend
> - Matrix assemble time by using 8 core : 70,234sec
>   int start, end;
>   VecGetOwnershipRange( element_vec, , );
>   for (int i=start; i~~ element matrix calculation
>matassemblybegin
>matassemblyend
> 
> 
> As you see the state, the parallel case spent a lot of time than sequential 
> case..
> How can I speed up in this case?
> Can I get some keyword or examples for parallelizing assembly of matrix in 
> finite element analysis ?
> 
> Thanks,
> Hyung Kim
> 



Re: [petsc-users] Union of sequential vecs

2022-12-12 Thread Karthikeyan Chockalingam - STFC UKRI via petsc-users
Thank you Matt and Blaise.

I will try out IS (though I have not it before)

(i)  Can IS be of different size, on different processors, and 
still call ISALLGather?

(ii)Can IS be passed as row indices to MatZeroRowsColumns?

I will look into DMlabels and start a different thread if needed.

Best,
Karthik.



From: Matthew Knepley 
Date: Saturday, 10 December 2022 at 02:20
To: Chockalingam, Karthikeyan (STFC,DL,HC) 
Cc: Mark Adams , petsc-users@mcs.anl.gov 

Subject: Re: [petsc-users] Union of sequential vecs
On Fri, Dec 9, 2022 at 6:50 PM Karthikeyan Chockalingam - STFC UKRI via 
petsc-users mailto:petsc-users@mcs.anl.gov>> wrote:
Thank you Mark and Barry.

@Mark Adams I follow you for the most part. Shouldn’t R 
be an MPI Vector?

Here it goes:

Q[Vec[i]] for all "i" in my local (sequential) "Vec".
Compute number of local nonzeros in Q, call it n.

//Create a MPI Vector R, with local size n


VecCreate(PETSC_COMM_WORLD, );

VecSetType(R, VECMPI);

VecSetSizes(R, n, PETSC_DECIDE);


//Populate MPI Vector R with local (sequential) "Vec".


VecGetOwnershipRange(R, , );

local_size = iend - istart; //The local_size should be ‘n’ right?

VecGetArray(R, );

  for (i = 0; i < local_size; i++) {

values[i] = Vec[i];

  }

VecRestoreArray(R, );


//Scatter R to all processors


VecV_SEQ;
VecScatter ctx;

VecScatterCreateToAll(R,,_SEQ);

//Remove duplicates in V_SEQ
How can I use PetscSortRemoveDupsInt to remove duplicates in V_SEQ?


Physics behind:
I am reading a parallel mesh, and want to mark all the boundary nodes. I use a 
local (sequential) Vec to store the boundary nodes for each parallel partition. 
Hence, local Vecs can end up with duplicate node index among them, which I 
would like to get rid of when I combine all of them together.

1) Blaise is right you should use an IS, not a Vec, to hold node indices. His 
solution is only a few lines, so I would use it.

2) I would not recommend doing things this way in the first place. PETSc can 
manage parallel meshes scalably, marking boundaries using DMLabel objects.

  Thanks,

 Matt

Best,
Karthik.





From: Mark Adams mailto:mfad...@lbl.gov>>
Date: Friday, 9 December 2022 at 21:08
To: Chockalingam, Karthikeyan (STFC,DL,HC) 
mailto:karthikeyan.chockalin...@stfc.ac.uk>>
Cc: Barry Smith mailto:bsm...@petsc.dev>>, 
petsc-users@mcs.anl.gov 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Union of sequential vecs
If your space is pretty compact, eg, (0,12), you could create an MPI vector Q 
of size 13, say, and each processor can add 1.0 to Q[Vec[i]], for all "i" in my 
local "Vec".
Then each processor can count the number of local nonzeros in Q, call it n, 
create a new vector, R, with local size n, then set R[i] = global index of the 
nonzero for each nonzero in Q, i=0:n.
Do some sort of vec-scatter-to-all with R to get what you want.

Does that work?

Mark


On Fri, Dec 9, 2022 at 3:25 PM Karthikeyan Chockalingam - STFC UKRI via 
petsc-users mailto:petsc-users@mcs.anl.gov>> wrote:
That is where I am stuck, I don’t know who to combine them to get Vec = 
{2,5,7,8,10,11,12}.
I just want them in an MPI vector.

I finally plan to call VecScatterCreateToAll so that all processor gets a copy.

Thank you.

Kind regards,
Karthik.

From: Barry Smith mailto:bsm...@petsc.dev>>
Date: Friday, 9 December 2022 at 20:04
To: Chockalingam, Karthikeyan (STFC,DL,HC) 
mailto:karthikeyan.chockalin...@stfc.ac.uk>>
Cc: petsc-users@mcs.anl.gov 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Union of sequential vecs

  How are you combining them to get Vec = {2,5,7,8,10,11,12}?

  Do you want the values to remain on the same MPI rank as before, just in an 
MPI vector?



On Dec 9, 2022, at 2:28 PM, Karthikeyan Chockalingam - STFC UKRI via 
petsc-users mailto:petsc-users@mcs.anl.gov>> wrote:

Hi,

I want to take the union of a set of sequential vectors, each living in a 
different processor.

Say,
Vec_Seq1 = {2,5,7}
Vec_Seq2 = {5,8,10,11}
Vec_Seq3 = {5,2,12}.

Finally, get the union of all them Vec = {2,5,7,8,10,11,12}.

I initially wanted to create a parallel vector and insert the (sequential 
vector) values but I do not know, to which index to insert the values to. But I 
do know the total size of Vec (which in this case is 7).

Any help is much appreciated.

Kind regards,
Karthik.




This email and any attachments are intended solely for the use of the named 
recipients. If you are not the intended recipient you must not use, disclose, 
copy or distribute this email or any of its attachments and should notify the 
sender immediately and delete this email from your system. UK Research and 
Innovation (UKRI) has taken every reasonable precaution to minimise risk of 
this email or any attachments containing viruses or malware but the recipient 
should carry out its own virus and malware checks before opening the 

Re: [petsc-users] parallelize matrix assembly process

2022-12-12 Thread 김성익
Hello Mark,

Following your comments,
I did run with '-info' and the outputs are as below
[image: image.png]
Global matrix seem to have preallocated well enough
And, ass I said earlier in the former email, If I run this code with mpi ,
It will be 70,000secs..
In this case, What is the problem?

And my loops already iterates over elements.
element_vec is just 1~125,000 array. for getting proper element index in
each process.
That example code is just simple schematic of my code.

Thanks,
Hyung Kim

2022년 12월 12일 (월) 오후 10:24, Mark Adams 님이 작성:

> Hi Hyung,
>
> First, verify that you are preallocating correctly.
> Run with '-info' and grep on "alloc" in the large output that you get.
> You will see lines like "number of mallocs in assembly: 0". You want 0.
> Do this with one processor and the 8.
>
> I don't understand your loop. You are iterating over vertices. You want to
> iterate over elements.
>
> Mark
>
>
>
> On Mon, Dec 12, 2022 at 6:16 AM 김성익  wrote:
>
>> Hello,
>>
>>
>> I need some keyword or some examples for parallelizing matrix assemble
>> process.
>>
>> My current state is as below.
>> - Finite element analysis code for Structural mechanics.
>> - problem size : 3D solid hexa element (number of elements : 125,000),
>> number of degree of freedom : 397,953
>> - Matrix type : seqaij, matrix set preallocation by using
>> MatSeqAIJSetPreallocation
>> - Matrix assemble time by using 1 core : 120 sec
>>for (int i=0; i<125000; i++) {
>> ~~ element matrix calculation}
>>matassemblybegin
>>matassemblyend
>> - Matrix assemble time by using 8 core : 70,234sec
>>   int start, end;
>>   VecGetOwnershipRange( element_vec, , );
>>   for (int i=start; i>~~ element matrix calculation
>>matassemblybegin
>>matassemblyend
>>
>>
>> As you see the state, the parallel case spent a lot of time than
>> sequential case..
>> How can I speed up in this case?
>> Can I get some keyword or examples for parallelizing assembly of matrix
>> in finite element analysis ?
>>
>> Thanks,
>> Hyung Kim
>>
>>


Re: [petsc-users] parallelize matrix assembly process

2022-12-12 Thread Mark Adams
Hi Hyung,

First, verify that you are preallocating correctly.
Run with '-info' and grep on "alloc" in the large output that you get.
You will see lines like "number of mallocs in assembly: 0". You want 0.
Do this with one processor and the 8.

I don't understand your loop. You are iterating over vertices. You want to
iterate over elements.

Mark



On Mon, Dec 12, 2022 at 6:16 AM 김성익  wrote:

> Hello,
>
>
> I need some keyword or some examples for parallelizing matrix assemble
> process.
>
> My current state is as below.
> - Finite element analysis code for Structural mechanics.
> - problem size : 3D solid hexa element (number of elements : 125,000),
> number of degree of freedom : 397,953
> - Matrix type : seqaij, matrix set preallocation by using
> MatSeqAIJSetPreallocation
> - Matrix assemble time by using 1 core : 120 sec
>for (int i=0; i<125000; i++) {
> ~~ element matrix calculation}
>matassemblybegin
>matassemblyend
> - Matrix assemble time by using 8 core : 70,234sec
>   int start, end;
>   VecGetOwnershipRange( element_vec, , );
>   for (int i=start; i~~ element matrix calculation
>matassemblybegin
>matassemblyend
>
>
> As you see the state, the parallel case spent a lot of time than
> sequential case..
> How can I speed up in this case?
> Can I get some keyword or examples for parallelizing assembly of matrix in
> finite element analysis ?
>
> Thanks,
> Hyung Kim
>
>


[petsc-users] parallelize matrix assembly process

2022-12-12 Thread 김성익
Hello,


I need some keyword or some examples for parallelizing matrix assemble
process.

My current state is as below.
- Finite element analysis code for Structural mechanics.
- problem size : 3D solid hexa element (number of elements : 125,000),
number of degree of freedom : 397,953
- Matrix type : seqaij, matrix set preallocation by using
MatSeqAIJSetPreallocation
- Matrix assemble time by using 1 core : 120 sec
   for (int i=0; i<125000; i++) {
~~ element matrix calculation}
   matassemblybegin
   matassemblyend
- Matrix assemble time by using 8 core : 70,234sec
  int start, end;
  VecGetOwnershipRange( element_vec, , );
  for (int i=start; i