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) >> >> >> > > > > -- 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)
