[petsc-users] Set matrix column to vector

[Florian Lindner]

Am 01.09.2014 12:45, schrieb Matthew Knepley:
> On Mon, Sep 1, 2014 at 4:10 AM, Florian Lindner <mailinglists at xgm.de> 
> wrote:
> 
>> Hello,
>> 
>> I want to set the entire column of a N x M matrix to a N vector. What 
>> is
>> the best way to do that?
>> 
>> My first guess would be to VecGetArray and use that array for
>> MatSetValuesLocal with nrow = VecGetLocalSize. What is the best to say
>> MatSetValuesLocal that I want to set all rows continuesly (same like
>> passing irow = [0, 1, ..., VecGetLocalSize-1]?
>> 
>> Any better way?
>> 
> 
> You are assuming dense storage above, so you can use
> 
> 
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatDenseGetArray.html+

How can you tell that I'm assuming dense storage. My matrix is actually 
dense, but I try to write my code as generic as possible (being a very 
petsc newbie). I have that code which crashes at the moment:

     void set_column_vector(Vector v, int col)
     {
       PetscErrorCode ierr = 0;
       const PetscScalar *vec;
       PetscInt size, mat_rows, mat_cols;
       VecGetLocalSize(v.vector, &size);
       cout << "Vector Size = " << size << endl;

       MatGetSize(matrix, &mat_rows, &mat_cols);
       cout << "Matrix Rows  = " << mat_rows << " Columns = " << mat_cols 
<< endl;
       PetscInt irow[size];
       for (int i = 0; i < size; i++) {
         irow[i] = i;
       }

       ierr = VecGetArrayRead(v.vector, &vec); CHKERRV(ierr);
       ierr = MatSetValuesLocal(matrix, size-1, irow, 1, &col, vec, 
INSERT_VALUES); CHKERRV(ierr);
       ierr = VecRestoreArrayRead(v.vector, &vec); CHKERRV(ierr);
       ierr = MatAssemblyBegin(matrix, MAT_FINAL_ASSEMBLY); 
CHKERRV(ierr);
       ierr = MatAssemblyEnd(matrix, MAT_FINAL_ASSEMBLY); CHKERRV(ierr);
     }

v.vector is a Vec, matrix is a Mat. col = 1. It's compiled with mpic++, 
but started without, just ./a.out. Output is:

Vector Size = 20
Matrix Rows  = 20 Columns = 5
[0]PETSC ERROR: 
------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or 
-on_error_attach_debugger
[0]PETSC ERROR: or see 
http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSC 
ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to 
find memory corruption errors
[0]PETSC ERROR: likely location of problem given in stack below
[0]PETSC ERROR: ---------------------  Stack Frames 
------------------------------------
[0]PETSC ERROR: Note: The EXACT line numbers in the stack are not 
available,
[0]PETSC ERROR:       INSTEAD the line number of the start of the 
function
[0]PETSC ERROR:       is given.
[0]PETSC ERROR: [0] MatSetValuesLocal line 1950 
/home/florian/software/petsc/src/mat/interface/matrix.c
[0]PETSC ERROR: --------------------- Error Message 
--------------------------------------------------------------
[0]PETSC ERROR: Signal received
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html 
for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.5.1, unknown
[0]PETSC ERROR: ./a.out on a arch-linux2-c-debug named asaru by florian 
Mon Sep  1 15:37:32 2014
[0]PETSC ERROR: Configure options --with-c2html=0
[0]PETSC ERROR: #1 User provided function() line 0 in  unknown file
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 59.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------


Any idea what the problem is?

Thanks!
Florian