You can use: MatSetOption(mat,MAT_IGNORE_ZERO_ENTRIES)
So subsequent MatSetValues() will ignore these zero entries. Satish On Wed, 14 Nov 2007, Tim Stitt wrote: > Matthew, > > OK...I see what you are saying. > > I initially set A a row at a time but for performance reasons I thought doing > it at once would be better. I overlooked the fact that the logical 2D matrix > input to MatSetValues() is non-zero values only. With hindsight I now remember > that was the case for each individual row. > > Thanks for pointing that out.... > > Regards. > > Matthew Knepley wrote: > > You appear to be setting every value in the sparse matrix. We do not > > throw out 0 values (since sometimes they are necessary for structural > > reasons). Thus you are allocating a ton of times. You need to remove > > the 0 values before calling MatSetValues (and their associated > > column entires as well). > > > > Matt > > > > On Nov 14, 2007 8:13 AM, Tim Stitt <timothy.stitt at ichec.ie> wrote: > > > > > Dear PETSc Users/Developers, > > > > > > I have the following sequential Fortran PETSc code that I have been > > > developing (on and off) based on the kind advice given by members of > > > this list, with respect to solving an inverse sparse matrix problem. > > > Essentially, the code reads in a square double complex matrix from > > > external file of size (order x order) and then proceeds to do a > > > MatMatSolve(), where A is the sparse matrix to invert, B is a dense > > > identity matrix and X is the resultant dense matrix....hope that makes > > > sense. > > > > > > My main problem is that the code stalls on the MatSetValues() for the > > > sparse matrix A. With a trivial test matrix of (224 x 224) the program > > > terminates successfully (by successfully I mean all instructions > > > execute...I am not interested in the validity of X right now). > > > Unfortunately, when I move up to a (2352 x 2352) matrix the > > > MatSetValues() routine for matrix A is still in progress after 15 > > > minutes on one processor of our AMD Opteron IBM Cluster. I know that > > > people will be screaming "preallocation"...but I have tried to take this > > > into account by running a loop over the rows in A and counting the > > > non-zero values explicitly prior to creation. I then pass this vector > > > into the creation routine for the nnz argument. For the large (2352 x > > > 2352) problem that seems to be taking forever to set...at most there are > > > only 200 elements per row that are non-zero according to the counts. > > > > > > Can anyone explain why the MatSetValues() routine is taking such a long > > > time. Maybe this expected for this specific task...although it seems > > > very long? > > > > > > I did notice that on the trivial (224 x 224) run that I was still > > > getting mallocs (approx 2000) for the A assembly when I used the -info > > > command line parameter. I thought that it should be 0 if my > > > preallocation counts were exact? Does this hint that I am doing > > > something wrong. I have checked the code but don't see any obvious > > > problems in the logic...not that means anything. > > > > > > I would be grateful if someone could advise on this matter. Also, if you > > > have a few seconds to spare I would be grateful if some experts could > > > scan the remaining logic of the code (not in fine detail) to make sure > > > that I am doing all that I need to do to get this calculation > > > working...assuming I can resolve the MatSetValues() problem. > > > > > > Once again many thanks in advance, > > > > > > Tim. > > > > > > ! Initialise the PETSc MPI Harness > > > call PetscInitialize(PETSC_NULL_CHARACTER,error);CHKERRQ(error) > > > > > > call MPI_COMM_SIZE(PETSC_COMM_SELF,processes,error);CHKERRQ(error) > > > call MPI_COMM_RANK(PETSC_COMM_SELF,ID,error);CHKERRQ(error) > > > > > > ! Read in Matrix > > > open(321,file='Hamiltonian.bin',form='unformatted') > > > read(321) order > > > if (ID==0) then > > > print * > > > print *,processes," Processing Elements being used" > > > print * > > > print *,"Matrix has order ",order," rows by ",order," columns" > > > print * > > > end if > > > > > > allocate(matrix(order,order)) > > > read(321) matrix > > > close(321) > > > > > > ! Allocate array for nnz > > > allocate(numberZero(order)) > > > > > > ! Count number of non-zero elements in each matrix row > > > do row=1,order > > > count=0 > > > do column=1,order > > > if (matrix(row,column).ne.(0,0)) count=count+1 > > > end do > > > numberZero(row)=count > > > end do > > > > > > ! Declare a PETSc Matrices > > > > > > call > > > MatCreateSeqAIJ(PETSC_COMM_SELF,order,order,PETSC_NULL_INTEGER,numberZero,A,error);CHKERRQ(error) > > > call > > > MatCreateSeqAIJ(PETSC_COMM_SELF,order,order,0,PETSC_NULL_INTEGER,factorMat,error);CHKERRQ(error) > > > call > > > MatCreateSeqDense(PETSC_COMM_SELF,order,order,PETSC_NULL_SCALAR,X,error);CHKERRQ(error) > > > call > > > MatCreateSeqDense(PETSC_COMM_SELF,order,order,PETSC_NULL_SCALAR,B,error);CHKERRQ(error) > > > > > > ! Set up zero-based array indexing for use in MatSetValues > > > allocate(columnIndices(order)) > > > > > > do column=1,order > > > columnIndices(column)=column-1 > > > end do > > > > > > ! Need to transpose values array as row-major arrays are used. > > > call > > > MatSetValues(A,order,columnIndices,order,columnIndices,transpose(matrix),INSERT_VALUES,error);CHKERRQ(error) > > > > > > ! Assemble Matrix A > > > call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error) > > > call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error) > > > > > > deallocate(matrix) > > > > > > ! Create Index Sets for Factorisation > > > call > > > ISCreateGeneral(PETSC_COMM_SELF,order,columnIndices,indexSet,error);CHKERRQ(error) > > > call MatFactorInfoInitialize(info,error);CHKERRQ(error) > > > call ISSetPermutation(indexSet,error);CHKERRQ(error) > > > call > > > MatLUFactorSymbolic(A,indexSet,indexSet,info,factorMat,error);CHKERRQ(error) > > > call MatLUFactorNumeric(A,info,factorMat,error);CHKERRQ(error) > > > > > > ! A no-longer needed > > > call MatDestroy(A,error);CHKERRQ(error) > > > > > > one=(1,0) > > > > > > ! Set Diagonal elements in Identity Matrix B > > > do row=0,order-1 > > > call MatSetValue(B,row,row,one,INSERT_VALUES,error);CHKERRQ(error) > > > end do > > > > > > ! Assemble B > > > call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error) > > > call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error) > > > > > > ! Assemble X > > > call MatAssemblyBegin(X,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error) > > > call MatAssemblyEnd(X,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error) > > > > > > ! Solve AX=B > > > call MatMatSolve(factorMat,B,X,error);CHKERRQ(error) > > > > > > ! Deallocate Storage > > > deallocate(columnIndices) > > > > > > call MatDestroy(factorMat,error);CHKERRQ(error) > > > call MatDestroy(B,error);CHKERRQ(error) > > > call MatDestroy(X,error);CHKERRQ(error) > > > > > > call PetscFinalize(error) > > > > > > -- > > > Dr. Timothy Stitt <timothy_dot_stitt_at_ichec.ie> > > > HPC Application Consultant - ICHEC (www.ichec.ie) > > > > > > Dublin Institute for Advanced Studies > > > 5 Merrion Square - Dublin 2 - Ireland > > > > > > +353-1-6621333 (tel) / +353-1-6621477 (fax) > > > > > > > > > > > > > > > > > > > >
