December 2010 - petsc-users - mailman3.cels.anl.gov

KSP solver and Distributed arrays
by Mohamad M. Nasr-Azadani 15 Dec '10

15 Dec '10

Hi guys, A simple question. Can I solve a linear system [A]{x} = {b} using KSP solvers using Matrix and rhs, solution vectors with different number of ghost nodes (width)? As an example, can it be possible to solve for vector {x} DA, width=3, STAR stencil, [A] DA, width=1, STAR stencil {b} DA, width=1, STAR stencil I am suspecting that it is not possible. Thanks, Mohamad

3 3

Updating the ghost nodes for distributed arrays
by Mohamad M. Nasr-Azadani 15 Dec '10

15 Dec '10

Hi guys, Is it possible to update the ghost values from a global to a local vector for distributed arrays when global and local vectors are not from the same DA, but the global vectors are the same? This is the the code that I have, (the only difference between the two DA's is the width. So, I am assuming that any global vector created based on those are going to be the same) G_data is created based on DA_3D, whereas L_data2 is created based on DA_3D2. Vec G_data, L_data; Vec G_data2, L_data2; ierr = DACreate3d(PCW, DA_NONPERIODIC, DA_STENCIL_STAR, NX, NY, NZ, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, width, PETSC_NULL, PETSC_NULL, PETSC_NULL, &DA_3D); ierr = DACreate3d(PCW, DA_NONPERIODIC, DA_STENCIL_STAR, NX, NY, NZ, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, width+2, PETSC_NULL, PETSC_NULL, PETSC_NULL, &DA_3D2); ierr = DACreateGlobalVector(DA_3D, &G_data); CHKERRQ(ierr); ierr = DACreateLocalVector(DA_3D, &L_data); CHKERRQ(ierr); ierr = DACreateGlobalVector(DA_3D2, &G_data2); CHKERRQ(ierr); ierr = DACreateLocalVector(DA_3D2, &L_data2); CHKERRQ(ierr); /* =====> Is this possible? */ ierr = DAGlobalToLocalBegin(DA_3D2, G_data, INSERT_VALUES, L_data2);CHKERRQ(ierr); ierr = DAGlobalToLocalEnd(DA_3D2, G_data, INSERT_VALUES, L_data2);CHKERRQ(ierr); Thanks, Mohamad

2 2

IS from DA by coordinates
by Ethan Coon 15 Dec '10

15 Dec '10

Hi all, Is there a cleaner way to create an IS to a global vector on a DA for a subset of nodes using a coordinate value than the following? (Pardon my pseudo-code combo of python and c and imprecise arguments) // get the coordinates of the nodes DAGetCoordinateDA(da, dac) DAGetCoordinates(da, vecc) DAVecGetArray(dac, vecc, vecc_a) // generate a one-dof da with no ghosts and the same parallel // structure as the da DAGetOwnershipRanges(da, lx[], ly[], lz[]) DACreate3D(comm, M,N,P,len(lx),len(ly),len(lz),1,0, lx[], ly[], \ lz[], da_one) // get the global indices of the one-dof da, noting that because // we set the stencil size to zero, the local array of global // indices is the same size as the local portion of the global array // of coordinates DAGetCorners(xs, ys, zs, xl, yl, zl) DAGetGlobalIndices(da_one, xl*yl*zl, indices[]) // loop over the local array returned by // DAGetGlobalIndices and the coordinates, // comparing to the test global_block_indices = [] for i in range(xs, xs+xl): for j in range(ys, ys+yl): for k in range(zs, zs+zl): if vecc_a[k,j,i,:] == whatever_coordinate: global_block_indices.append(indices[k,j,i]) // make the IS ISCreateBlock(comm, ndofs, global_block_indices, coord_is) // restore and destroy etc. This just seems quite complicated with the construction of the one-dof da to get global indices of the block. There might be a better way with ISLocalToGlobalMapping, but I wasn't sure how. Any suggestions? Thanks, Ethan -- ------------------------------------- Ethan Coon Post-Doctoral Researcher Mathematical Modeling and Analysis Los Alamos National Laboratory 505-665-8289 http://www.ldeo.columbia.edu/~ecoon/ -------------------------------------

2 2

MatMult
by Jakub Pola 13 Dec '10

13 Dec '10

Hello again, I compiled one of te examples. I used sparse matix called 02-raefsky3. I used -vec_type cuda and -mat_type seqaijcuda. When I see summary of the operations performed by program there is MatMult 1 1.0 2.0237e-02 1.0 2.98e+06 1.0 0.0e+00 0.0e+00 0.0e+00 2100 0 0 0 2100 0 0 0 147 Does time of performing MatMult includes memory transfer for loading matrix in GPU memory or just exact computation time? Thanks in advance. Kuba.

3 6

class templates in C++ and Petsc functions
by Umut Tabak 13 Dec '10

13 Dec '10

Dear all, I was trying to write a simple class template code to interface Petsc matrices along with Boost sparse matrices and vice versa. However, my first naive try did not work out of the box, since the template class requires some functions from PETSc libraries and these functions should be located in the object files to be able to compile. Since class template is header only and there is no object code generation, it can not find the PETSc library functions which is logical. I am also new to templates in C++. I just gave a try to create a template class instead of code duplication however it did not end up as expected. Are there some smarter ways to accomplish this task. The function analyzes a boost csr matrix and sets the rows of a Petsc matrix. T is for boost matrices and T1 for Mat in Petsc. I can also get some comments on the code ;) Best regards, Umut template <class T, class T1> int converter<T, T1>::convertMBo2Pe( const T& A, T1 A_ ){ PetscErrorCode ierr; int cntNnz = 0; typedef typename T::iterator1 i1_t; typedef typename T::iterator2 i2_t; //int nnz[ooFelieMatrix.size1()]; int nnz[A.size1()]; unsigned ind=0; //get information about the matrix double* vals = NULL; for (i1_t i1 = A.begin1(); i1 != A.end1(); ++i1) { nnz[ind] = distance(i1.begin(), i1.end()); ind++; } // create the matrix depending // on the values of the nonzeros // on each row ierr = MatCreateSeqAIJ( PETSC_COMM_SELF, A.size1(), A.size2(), cntNnz, nnz, A_ ); PetscInt rInd = 0, cInd=0; PetscInt* rCount, dummy; rCount = &dummy; // pointer to values in a row PetscScalar* valsOfRowI = NULL; PetscInt* colIndexOfRowI = NULL; PetscInt rC = 1; for(i1_t i1 = A.begin1(); i1 != A.end1(); ++i1) { // allocate space for the values of row I valsOfRowI = new PetscScalar[nnz[rInd]]; colIndexOfRowI = new PetscInt[nnz[rInd]]; for(i2_t i2 = i1.begin(); i2 != i1.end(); ++i2) { colIndexOfRowI[cInd] = i2.index2(); valsOfRowI[cInd] = *i2; cInd++; } // setting one row each time *rCount = rInd; MatSetValues( A_, rC, rCount, nnz[rInd], colIndexOfRowI, valsOfRowI, INSERT_VALUES ); // delete delete [] valsOfRowI; delete [] colIndexOfRowI; rInd++; cInd = 0; } // MatAssemblyBegin( A_, MAT_FINAL_ASSEMBLY ); MatAssemblyEnd( A_, MAT_FINAL_ASSEMBLY ); // return return 0; } -- - Hope is a good thing, maybe the best of things and no good thing ever dies... The Shawshank Redemption, replique of Tim Robbins

2 1

Solution output of TS
by Michael E Henderson 13 Dec '10

13 Dec '10

Hi, I'm using TS and seeing a ":flickering" in the output. I use ierr=TSCreate(PETSC_COMM_WORLD,&timeStepper); ierr=TSSetProblemType(timeStepper,TS_NONLINEAR); ierr=TSSetFromOptions(timeStepper); ierr=TSSetIFunction(timeStepper, formOperatorImplicitTimeFunction , (void*)Ldata); ierr=TSSetIJacobian(timeStepper,A,A, formOperatorImplicitTimeJacobian , (void*)Ldata); ierr=TSMonitorSet(timeStepper,SaveTSSolution,(void*)(&data),NULL); I'm using the implicit nonlinear formulation for an ADE, and writing the solution in the monitoring routine. If I instead write out the solution in the IFunction routine when the residual is small I see a different (better) solution passed in. Does this sound like a problem that's been seen before? Thanks, MIke Henderson ------------------------------------------------------------------------------------------------------------------------------------ Mathematical Sciences, TJ Watson Research Center mhender(a)watson.ibm.com http://www.research.ibm.com/people/h/henderson/ http://multifario.sourceforge.net/

1 0

Writing PETSc matrices
by Nunion 13 Dec '10

13 Dec '10

Hello, I am new to PETSc and programming. I have a question concerning writing PETSc matrices in binary from binary matrices [compressessed/uncompressed] generated in Matlab. I am attempting to use the files in the /bin/matlab directory, in particular the PetscBinaryWrite.m file. However, the usage; PetscBinaryWrite('matrix.mat','output.ex') does not seem to work. I also tried using the examples in the /mat directory however, matlab does not support the writing of complex matrices in ASCII. Thanks in advance, Tom

5 6

Re: [petsc-users] run direct linear solver in parallel
by Max Ng 13 Dec '10

13 Dec '10

Hi, I am having a similar problem and I'm using PETSc 3.1-p6. I wish to use SPOOLES because I need to build on Windows with VC++ (and without a Fortran compiler). And in my tests somehow SPOOLES performs better than SuperLU. My program runs correctly in mpiexec -n 1. When I try mpiexec -n 2, I got this error: Assertion failed in file helper_fns.c at line 337: 0 memcpy argument memory ranges overlap, dst_=0x972ef84 src_=0x972ef84 len_=4 internal ABORT - process 1 Assertion failed in file helper_fns.c at line 337: 0 memcpy argument memory ranges overlap, dst_=0x90c4018 src_=0x90c4018 len_=4 internal ABORT - process 0 rank 1 in job 113 vm1_57881 caused collective abort of all ranks exit status of rank 1: killed by signal 9 Here is the source code: // N = 40000, n = 20000, nnz = 9 // MatCreate(comm, &mat); MatSetType(mat, MATAIJ); MatSetSizes(mat, n, n, N, N); MatSeqAIJSetPreallocation(mat, nnz, PETSC_NULL); MatMPIAIJSetPreallocation(mat, nnz, PETSC_NULL, nnz, PETSC_NULL); // some code to fill the matrix values // ... KSPCreate(comm, &ksp); KSPSetOperators(ksp, mat, mat, DIFFERENT_NONZERO_PATTERN); KSPSetType(ksp, KSPPREONLY); KSPGetPC(ksp, &pc); PCSetType(pc, PCLU); PCFactorSetMatSolverPackage(pc, MAT_SOLVER_SPOOLES); KSPSetUp(ksp); It crashes at the KSPSetUp() statement. Do you have any ideas? Thanks in advance! Max Ng On Dec 3, 2010, at 4:19 PM, Xiangdong Liang wrote: > Hi everyone,> > I am wondering how I can run the direct solver in parallel. I can run> my program in a single processor with direct linear solver by> > ./foo.out -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package spooles> > However, when I try to run it with mpi:> > mpirun.openmpi -np 2 ./foo.out -ksp_type preonly -pc_type lu> -pc_factor_mat_solver_package spooles> > I got error like this:> > [0]PETSC ERROR: --------------------- Error Message> ------------------------------------> [0]PETSC ERROR: No support for this operation for this object type!> [0]PETSC ERROR: Matrix type mpiaij symbolic LU!> > [0]PETSC ERROR: Libraries linked from> /home/hazelsct/petsc-2.3.3/lib/linux-gnu-c-opt> [0]PETSC ERROR: Configure run at Mon Jun 30 14:37:52 2008> [0]PETSC ERROR: Configure options --with-shared --with-dynamic> --with-debugging=0 --useThreads 0 --with-mpi-dir=/usr/lib/openmpi> --with-mpi-shared=1 --with-blas-lib=-lblas --with-lapack-lib=-llapack> --with-umfpack=1 --with-umfpack-include=/usr/include/suitesparse> --with-umfpack-lib="[/usr/lib/libumfpack.so,/usr/lib/libamd.so]"> --with-superlu=1 --with-superlu-include=/usr/include/superlu> --with-superlu-lib=/usr/lib/libsuperlu.so --with-spooles=1> --with-spooles-include=/usr/include/spooles> --with-spooles-lib=/usr/lib/libspooles.so --with-hypre=1> --with-hypre-dir=/usr --with-babel=1 --with-babel-dir=/usr> [0]PETSC ERROR:> ------------------------------------------------------------------------> [0]PETSC ERROR: MatLUFactorSymbolic() line 2174 in src/mat/interface/matrix.c> [0]PETSC ERROR: PCSetUp_LU() line 257 in src/ksp/pc/impls/factor/lu/lu.c> -------------------------------------------------------> > Would you like to tell me where I am doing wrong? I appreciate your help.> > Xiangdong

3 3

About the "PetscOptionsSetValue" usage
by Filippo Spiga 13 Dec '10

13 Dec '10

(sorry in advance for the cross-posting) Dear all, I recently decided to use in my program the routine "PetscOptionsSetValue" to allow to change at runtime the parameter of KSP/SNES. Of course, my code performs several operations on different KSP and SNES objects. Every object should have its own set of options (different tollerancies, different preconditioner) and I want to tune it at runtime. So "PetscOptionsSetValue" it is great. But I have a question: when I set an option, this option is still valid for all the program. Let's assume this situation (it is an example, I use fake options): PetscOptionsSetValue("A", "1"); PetscOptionsSetValue("B", "2"); PetscOptionsSetValue("C", "3"); ierr = KSPSolve(); PetscOptionsSetValue("A", "10"); PetscOptionsSetValue("B", "20"); ierr = SNESSolve(); This is what happens (if I correcly understood): - KSP runs with this presets {A=1, B=2, C=3} that are different than the defaults. - SNES runs with this presets {A=10, B=20, C=3}. So SNES runs with C=3 that is different from the default. But I would like to use the default because C=3 produces wrong errors. How I can easily reset the options' database? Thanks in advance, Cheers -- Filippo SPIGA, MSc Computer Science «Nobody will drive us out of Cantor's paradise.» -- David Hilbert ***** Disclaimer: "Please note this message and any attachment are CONFIDENTIAL an may be privileged or otherwise protected from disclosure. The contents are not to be disclosed to anyone other than the addressee. Unauthorized recipients are requested to preserve this confidentiality and to advise the sender immediately of any error in transmission."

2 5

Create csr matrix
by Jakub Pola 12 Dec '10

12 Dec '10

Hi, Could you please help me with creating CSR matrix. I have only one processor so It have to be done locally Here is the matrix from those information I would like to have matrix petsc matrix. double vals [] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0} ; int c_idx [] = {1, 3, 2, 0, 1, 2, 2, 3} ; int r_idx [] = {0, 2, 3, 6, 8} ; int n_rows = 4 ; //square matrix Thank you in advance for help. Kuba

2 4