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

potential e-mail disruption due to power outage
by Satish Balay 03 Dec '10

03 Dec '10

We are having a power outage this weekend at MCS [staring friday evening, and ending sunday evening or monday morning CST] e-mail, mailing-lists, ftp, web servers are supposed to work during this outage. But if something isn't working - and we are unable to respond by e-mail, you now know the reason. Satish

1 0

Shell matrix and MatMultAdd.
by Vijay S. Mahadevan 02 Dec '10

02 Dec '10

Hi all, There is probably some minor inconsistency in my understanding but is there a fundamental difference between the following two options ? // Option 1 ierr = MatMult(A, solution, temporaryvec) ;CHKERRQ(ierr); ierr = VecAXPY(rhs, 1.0, temporaryvec) ;CHKERRQ(ierr); // Option 2 ierr = MatMultAdd(A, solution, rhs, rhs) ;CHKERRQ(ierr); Here A is a shell matrix (serial) that has a routine defined to perform MATOP_MULT only. I ask because Option 1 gives me the right result while Option 2 segfaults at the MatMultAdd line. My only logical conclusion is that I need to define a MATOP_MULT_ADD or something similar for the shell for this to work. Is this understanding correct ? I implicitly assumed that petsc recognizes that MATOP_MULT has been defined already and since MatMultAdd only requires the action of a matrix on a vector to perform its operation, should this not be computed by petsc automatically ? I really do not want to creat an extra vector here with Option 1 since this occurs at a finer level in my calculation. But is there any way that you would suggest I do this without extra allocations ? Any comments or pointers will be much appreciated. Vijay

2 4

column index in MatSetValues()
by Peter Wang 02 Dec '10

02 Dec '10

I am trying to create a matrix and insert values to it. The martix is supposed to be as following: 1 0 0 0 0 2 0 0 0 0 3 0 0 0 0 4 array coef[] is the diagonal value of the matrix, snr[] is the index of the row, rnr[] is the index of column. However, I always get the wrong results. It shows the Column too large: col 4607182418800017408 max 3! I cheked the value of rnr[]. The output snr and rnr is correct: snr= 0 1 2 3 rnr= 0 1 2 3 It seems there is something wrong when MatSetValues() is called. Following is a part of the error information. The information is shown at each loop of do II=Istart,Iend-1 The output (if any) follows: snr= 0 1 2 3 rnr= 0 1 2 3 8.....Check after MatGetOwnershipRange() Istart= 0 Iend= 4 II= 0 1 0 0 [0]PETSC ERROR: --------------------- Error Message ------------------------------------ [0]PETSC ERROR: Argument out of range! [0]PETSC ERROR: Column too large: col 4607182418800017408 max 3! [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: Petsc Release Version 3.1.0, Patch 5, Mon Sep 27 11:51:54 CDT 2010 [0]PETSC ERROR: See docs/changes/index.html for recent updates. [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting. [0]PETSC ERROR: See docs/index.html for manual pages. [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: Debug_PETSc_MatCreate_20101130 on a linux-gnu named compute-1-35.hpc.local.uwm by pwang_a Tue Nov 30 22:27:03 2010 [0]PETSC ERROR: Libraries linked from /sharedapps/uwm/ceas/gcc-4.4.3/petsc/3.1-p5-v1/lib [0]PETSC ERROR: Configure run at Fri Oct 8 12:59:16 2010 [0]PETSC ERROR: Configure options --prefix=/sharedapps/uwm/ceas/gcc-4.4.3/petsc/3.1-p5-v1 --with-mpi-dir=/sharedapps/uwm/common/gcc-4.4.3/openmpi/1.3.2-v1 --with-blas-lapack-dir=/sharedapps/uwm/ceas/gcc-4.4.3/lapack/3.2.2-v1/lib --with-64-bit-indices=1 --with-64-bit-pointers=1 --with-large-file-io=1 --with-x=0 [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: MatSetValues_SeqAIJ() line 193 in src/mat/impls/aij/seq/aij.c [0]PETSC ERROR: MatSetValues() line 992 in src/mat/interface/matrix.c !The code is as following: !============================= program Debug_PETSc_MatCreate_20101130 implicit none ! #include "finclude/petscsys.h" #include "finclude/petscvec.h" #include "finclude/petscmat.h" #include "finclude/petscpc.h" #include "finclude/petscksp.h" ! Variables !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! PETSc Variables !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - real*8 norm PetscInt i,j,II,JJ,its !,m,n PetscInt Istart,Iend,ione PetscErrorCode ierr PetscMPIInt myid,numprocs PetscTruth flg PetscScalar v,one,neg_one Vec x,b,u Mat A_petsc KSP ksp PetscInt,parameter:: n_nz=4 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Other Variables !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - !parameter::n_nz=4 Real*8::Coef(n_nz) PetscInt::snr(n_nz),rnr(n_nz) data Coef /1., 2., 3. , 4./ data snr /0, 1, 2, 3/ data rnr /0, 1 , 2, 3/ ! Body of Debug_PETSc_MatCreate_20101130 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Beginning of program ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - call PetscInitialize(PETSC_NULL_CHARACTER,ierr) call MPI_Comm_rank(PETSC_COMM_WORLD,myid,ierr) call MPI_Comm_size(PETSC_COMM_WORLD,numprocs,ierr) write(*,"('snr=',4i4)")snr write(*,"('rnr=',4i4)")rnr call MatCreate(PETSC_COMM_WORLD,A_Petsc,ierr) call MatSetSizes(A_Petsc,PETSC_DECIDE,PETSC_DECIDE,n_nz,n_nz,ierr) !n_nz-1??? call MatSetFromOptions(A_Petsc,ierr) ! write(*,*)A_petsc call MatGetOwnershipRange(A_Petsc,Istart,Iend,ierr) write(*,'(1a,1i7,1a,1i7)') & '8.....Check after MatGetOwnershipRange() Istart=',Istart,' Iend=',Iend do II=Istart,Iend-1 ione=II+1 !(Coef,snr,rnr are 1-based row and column numbers, shifting them to 0-based) write(*,'(1a,4i7)')'II=',II,ione,snr(ione),rnr(ione) !output snr and rnr for error check call MatSetValues(A_Petsc,1,snr(Ione),1,rnr(Ione),Coef(Ione),INSERT_VALUES,ierr) enddo write(*,'(1a)')'9.....Check after MatSetValues()' call MatAssemblyBegin(A_petsc,MAT_FINAL_ASSEMBLY,ierr) call MatAssemblyEnd(A_Petsc,MAT_FINAL_ASSEMBLY,ierr) write(*,'(1a)')'10.....Check after MatCreate()' call MatView(A_Petsc,PETSC_VIEWER_STDOUT_WORLD,ierr) ! call KSPDestroy(ksp,ierr) ! call VecDestroy(u,ierr) ! call VecDestroy(x,ierr) ! call VecDestroy(b,ierr) call MatDestroy(A_petsc,ierr) call PetscFinalize(ierr) end program Debug_PETSc_MatCreate_20101130 !=====================================

2 7

What's the binary function corresponding to "PetscViewerASCIIPrintf"?
by Li, Zhisong (lizs) 01 Dec '10

01 Dec '10

Hi, Petsc Team, I once read that when we write large volume of data into an output data file, it's best to write in binary format rather than in ASCII format. I forget the origin of this statement and could not find it any more. What I found in the example code is about writing array into a binary file. Besides data stored in arrays, I also want to include some other info such as document title, data dimensions and line changing, which are not arrays, into the output file. So a "PetscViewerASCIIPrintf" is very convenient for doing this if we use ASCII format. But I could not find its correspondent function for binary format. Can you tell me how to do this? A sample code will be most helpful. Best Regards, Zhisong Li

2 1

Why MatMult(A, x, Ax) != A*x?
by Rongliang Chen 30 Nov '10

30 Nov '10

Hi, I have tried the following method: Vec Ax; ierr = VecDuplicate(x, &Ax);CHKERRQ(ierr); ierr = MatMult(*A, x, Ax);CHKERRQ(ierr); PetscViewer viewer = PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD); ierr = MatView(*A,viewer);CHKERRQ(ierr); ierr = VecView(x,viewer);CHKERRQ(ierr); ierr = VecView(Ax,viewer);CHKERRQ(ierr); // ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr); (If I use this line, I just can obtain A,x,y and I can not get the A1, x1, y1 (see below.)). After do two SNES, I use the following commands in matlab: >> [A,x,y,A1,x1,y1] = PetscBinaryRead('binaryoutput'); >> norm(A*x-y) ans = 1.2000e-13 >> norm(A1*x1-y1) ans = 2.6612e+03 The reason that I want to do this is that I doubt that there is a bug in my code and I want to following these information to find out the bug. There is a very strange thing in my code:(I use GMRES with ASM preconditioner to solve the linear system.) Solving Shape Optimization problem 0 SNES norm 1.1359606379e+02, 0 KSP its last norm 0.0000000000e+00. 1 SNES norm 4.1335342946e+01, 1 KSP its last norm 2.3855348597e-03. 2 SNES norm *4.1085343652e+01, 13 KSP* its last norm 3.4528582271e-02. *SNES diverged DIVERGED_LS_FAILURE.* I think there is some thing wrong with the second SNES, because its KSP should be 1 or 2. If I destroy the Jacobian matrix here and recreate it, then continue the SNES, it will be alright: 0 SNES norm *4.1085343652e+01,* 0 KSP its last norm 0.0000000000e+00. 1 SNES norm 6.9950501590e+00, 1 KSP its last norm 2.8278185399e-03. 2 SNES norm 1.3445253365e-01, 1 KSP its last norm 2.3243944231e-04. 3 SNES norm 5.0914645029e-04, 1 KSP its last norm 5.0954271956e-06. ................. So I output the matrix and vectors into matlab and do the matrix and vector multiplication and I found that the result is different from the MatMult in Petsc. Best, Rongliang ------------------------------ > > Message: 2 > Date: Tue, 30 Nov 2010 21:52:45 +0100 > From: Jed Brown <jed(a)59A2.org> > Subject: Re: [petsc-users] Why MatMult(A, x, Ax) != A*x? > To: PETSc users list <petsc-users(a)mcs.anl.gov> > Message-ID: > <AANLkTiks1Go4J56qYDuFabb+QQBtx0DmfaythUkzgJik(a)mail.gmail.com<AANLkTiks1Go4J56qYDuFabb%2BQQBtx0DmfaythUkzgJik(a)mail.gmail.com> > > > Content-Type: text/plain; charset="utf-8" > > On Tue, Nov 30, 2010 at 21:42, Rongliang Chen <rongliang.chan(a)gmail.com > >wrote: > > > I know that I should call PetscViewerDestroy() when I have done with each > > viewer. But if I called this I just can view the last Matrix and vector. > > What I want > > to do is that I let it do three steps of SNES, view the jacobian of each > > step. > > > > Maybe view it to different files. The code you posted for viewing is > certainly not correct, especially if called in a loop. > > > > I found that A*x is correct at the first step and it's wrong for the > second > > and > > third step. > > I also tried using the binary output and input, and it has the same > > results. > > > > The problem is most likely that you are picking these pieces out of the > files differently than they are put in. MatMult being incorrect is by far > the least likely scenario. I don't know why you want to do this, but you > need to find a way to manage the files so that you can be certain that you > are reading (from Matlab) exactly what you think you are reading. One way > to do this is to use a binary file that gets all three items together. > > Jed > > >

3 2