On Tue, Jul 3, 2012 at 6:18 AM, Klaij, Christiaan <span dir="ltr"><<a href="mailto:C.Klaij@marin.nl" target="_blank">C.Klaij@marin.nl</a>></span> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
I'm trying to understand the use of null spaces. Whatever I do, it<br>
always seem to pass the null space test. Could you please tell me<br>
what's wrong with this example (c++, petsc-3.3-p1):<br></blockquote><div><br></div><div>Yes, there was a bug in 3.3. I pushed a fix which should go out in the next patch, and</div><div>its in petsc-dev.</div><div><br>
</div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
$ cat nullsp.cc<br>
// test null space check<br>
<br>
#include <petscksp.h><br>
<br>
int main(int argc, char **argv) {<br>
<br>
PetscErrorCode ierr;<br>
PetscInt row,start,end;<br>
PetscScalar val[1];<br>
PetscReal norm;<br>
Mat A;<br>
Vec x,y;<br>
MatNullSpace nullsp;<br>
PetscBool isNull;<br>
<br>
ierr = PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL); CHKERRQ(ierr);<br>
<br>
// diagonal matrix<br>
ierr = MatCreate(PETSC_COMM_WORLD,&A); CHKERRQ(ierr);<br>
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,24,24); CHKERRQ(ierr);<br>
ierr = MatSetType(A,MATMPIAIJ); CHKERRQ(ierr);<br>
ierr = MatMPIAIJSetPreallocation(A,1,PETSC_NULL,1,PETSC_NULL); CHKERRQ(ierr);<br>
ierr = MatGetOwnershipRange(A,&start,&end); CHKERRQ(ierr);<br>
for (row=start; row<end; row++) {<br>
val[0] = 1.0;<br>
ierr = MatSetValues(A,1,&row,1,&row,val,INSERT_VALUES); CHKERRQ(ierr);<br>
}<br>
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
<br>
// random vector<br>
ierr = VecCreate(PETSC_COMM_WORLD,&x); CHKERRQ(ierr);<br>
ierr = VecSetSizes(x,PETSC_DECIDE,24); CHKERRQ(ierr);<br>
ierr = VecSetType(x,VECMPI); CHKERRQ(ierr);<br>
ierr = VecSetRandom(x,PETSC_NULL); CHKERRQ(ierr);<br>
<br>
// null space<br>
ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,1,&x,&nullsp); CHKERRQ(ierr);<br>
ierr = MatNullSpaceTest(nullsp,A,&isNull); CHKERRQ(ierr);<br>
if (isNull==PETSC_TRUE) {<br>
ierr = PetscPrintf(PETSC_COMM_WORLD,"null space check passed\n");<br>
}<br>
else {<br>
ierr = PetscPrintf(PETSC_COMM_WORLD,"null space check failed\n");<br>
}<br>
<br>
// check null space<br>
ierr = VecDuplicate(x,&y); CHKERRQ(ierr);<br>
ierr = MatMult(A,x,y); CHKERRQ(ierr);<br>
ierr = VecNorm(y,NORM_2,&norm); CHKERRQ(ierr);<br>
ierr = PetscPrintf(PETSC_COMM_WORLD,"|Ax| = %G\n",norm); CHKERRQ(ierr);<br>
<br>
ierr = PetscFinalize(); CHKERRQ(ierr);<br>
<br>
return 0;<br>
}<br>
<br>
$ mpiexec -n 2 ./nullsp -log_summary<br>
null space check passed<br>
|Ax| = 2.51613<br>
************************************************************************************************************************<br>
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***<br>
************************************************************************************************************************<br>
<br>
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------<br>
<br>
./nullsp on a linux_64b named lin0133 with 2 processors, by cklaij Tue Jul 3 14:03:34 2012<br>
Using Petsc Release Version 3.3.0, Patch 1, Fri Jun 15 09:30:49 CDT 2012<br>
<br>
Max Max/Min Avg Total<br>
Time (sec): 4.794e-03 1.00000 4.794e-03<br>
Objects: 1.400e+01 1.00000 1.400e+01<br>
Flops: 7.200e+01 1.00000 7.200e+01 1.440e+02<br>
Flops/sec: 1.502e+04 1.00000 1.502e+04 3.004e+04<br>
Memory: 5.936e+04 1.00000 1.187e+05<br>
MPI Messages: 0.000e+00 0.00000 0.000e+00 0.000e+00<br>
MPI Message Lengths: 0.000e+00 0.00000 0.000e+00 0.000e+00<br>
MPI Reductions: 4.100e+01 1.00000<br>
<br>
Flop counting convention: 1 flop = 1 real number operation of type (multiply/divide/add/subtract)<br>
e.g., VecAXPY() for real vectors of length N --> 2N flops<br>
and VecAXPY() for complex vectors of length N --> 8N flops<br>
<br>
Summary of Stages: ----- Time ------ ----- Flops ----- --- Messages --- -- Message Lengths -- -- Reductions --<br>
Avg %Total Avg %Total counts %Total Avg %Total counts %Total<br>
0: Main Stage: 4.7846e-03 99.8% 1.4400e+02 100.0% 0.000e+00 0.0% 0.000e+00 0.0% 4.000e+01 97.6%<br>
<br>
------------------------------------------------------------------------------------------------------------------------<br>
See the 'Profiling' chapter of the users' manual for details on interpreting output.<br>
Phase summary info:<br>
Count: number of times phase was executed<br>
Time and Flops: Max - maximum over all processors<br>
Ratio - ratio of maximum to minimum over all processors<br>
Mess: number of messages sent<br>
Avg. len: average message length<br>
Reduct: number of global reductions<br>
Global: entire computation<br>
Stage: stages of a computation. Set stages with PetscLogStagePush() and PetscLogStagePop().<br>
%T - percent time in this phase %f - percent flops in this phase<br>
%M - percent messages in this phase %L - percent message lengths in this phase<br>
%R - percent reductions in this phase<br>
Total Mflop/s: 10e-6 * (sum of flops over all processors)/(max time over all processors)<br>
------------------------------------------------------------------------------------------------------------------------<br>
<br>
<br>
##########################################################<br>
# #<br>
# WARNING!!! #<br>
# #<br>
# This code was compiled with a debugging option, #<br>
# To get timing results run ./configure #<br>
# using --with-debugging=no, the performance will #<br>
# be generally two or three times faster. #<br>
# #<br>
##########################################################<br>
<br>
<br>
Event Count Time (sec) Flops --- Global --- --- Stage --- Total<br>
Max Ratio Max Ratio Max Ratio Mess Avg len Reduct %T %f %M %L %R %T %f %M %L %R Mflop/s<br>
------------------------------------------------------------------------------------------------------------------------<br>
<br>
--- Event Stage 0: Main Stage<br>
<br>
MatMult 1 1.0 4.0054e-05 1.1 1.20e+01 1.0 0.0e+00 0.0e+00 0.0e+00 1 17 0 0 0 1 17 0 0 0 1<br>
MatAssemblyBegin 1 1.0 3.7909e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 2.0e+00 1 0 0 0 5 1 0 0 0 5 0<br>
MatAssemblyEnd 1 1.0 4.3201e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 1.9e+01 9 0 0 0 46 9 0 0 0 48 0<br>
VecNorm 2 1.0 2.9700e-03 1.0 4.80e+01 1.0 0.0e+00 0.0e+00 2.0e+00 62 67 0 0 5 62 67 0 0 5 0<br>
VecSet 1 1.0 5.0068e-06 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0<br>
VecScatterBegin 2 1.0 7.8678e-06 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0<br>
VecScatterEnd 2 1.0 4.0531e-06 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0<br>
VecSetRandom 1 1.0 9.0599e-06 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0<br>
------------------------------------------------------------------------------------------------------------------------<br>
<br>
Memory usage is given in bytes:<br>
<br>
Object Type Creations Destructions Memory Descendants' Mem.<br>
Reports information only for process 0.<br>
<br>
--- Event Stage 0: Main Stage<br>
<br>
Matrix 3 0 0 0<br>
Matrix Null Space 1 0 0 0<br>
Vector 5 1 1504 0<br>
Vector Scatter 1 0 0 0<br>
Index Set 2 2 1496 0<br>
PetscRandom 1 1 616 0<br>
Viewer 1 0 0 0<br>
========================================================================================================================<br>
Average time to get PetscTime(): 9.53674e-08<br>
Average time for MPI_Barrier(): 4.29153e-07<br>
Average time for zero size MPI_Send(): 8.58307e-06<br>
#PETSc Option Table entries:<br>
-log_summary<br>
#End of PETSc Option Table entries<br>
Compiled without FORTRAN kernels<br>
Compiled with full precision matrices (default)<br>
sizeof(short) 2 sizeof(int) 4 sizeof(long) 8 sizeof(void*) 8 sizeof(PetscScalar) 8 sizeof(PetscInt) 4<br>
Configure run at: Wed Jun 20 12:08:20 2012<br>
Configure options: --with-mpi-dir=/opt/refresco/libraries_cklaij/openmpi-1.4.5 --with-clanguage=c++ --with-x=1 --with-debugging=1 --with-hypre-include=/opt/refresco/libraries_cklaij/hypre-2.7.0b/include --with-hypre-lib=/opt/refresco/libraries_cklaij/hypre-2.7.0b/lib/libHYPRE.a --with-ml-include=/opt/refresco/libraries_cklaij/ml-6.2/include --with-ml-lib=/opt/refresco/libraries_cklaij/ml-6.2/lib/libml.a --with-blas-lapack-dir=/opt/intel/mkl<br>
-----------------------------------------<br>
Libraries compiled on Wed Jun 20 12:08:20 2012 on lin0133<br>
Machine characteristics: Linux-2.6.32-41-generic-x86_64-with-Ubuntu-10.04-lucid<br>
Using PETSc directory: /home/CKlaij/ReFRESCO/Dev/trunk/Libs/build/petsc-3.3-p1<br>
Using PETSc arch: linux_64bit_debug<br>
-----------------------------------------<br>
<br>
Using C compiler: /opt/refresco/libraries_cklaij/openmpi-1.4.5/bin/mpicxx -wd1572 -g ${COPTFLAGS} ${CFLAGS}<br>
Using Fortran compiler: /opt/refresco/libraries_cklaij/openmpi-1.4.5/bin/mpif90 -g ${FOPTFLAGS} ${FFLAGS}<br>
-----------------------------------------<br>
<br>
Using include paths: -I/home/CKlaij/ReFRESCO/Dev/trunk/Libs/build/petsc-3.3-p1/linux_64bit_debug/include -I/home/CKlaij/ReFRESCO/Dev/trunk/Libs/build/petsc-3.3-p1/include -I/home/CKlaij/ReFRESCO/Dev/trunk/Libs/build/petsc-3.3-p1/include -I/home/CKlaij/ReFRESCO/Dev/trunk/Libs/build/petsc-3.3-p1/linux_64bit_debug/include -I/opt/refresco/libraries_cklaij/hypre-2.7.0b/include -I/opt/refresco/libraries_cklaij/ml-6.2/include -I/opt/refresco/libraries_cklaij/openmpi-1.4.5/include<br>
-----------------------------------------<br>
<br>
Using C linker: /opt/refresco/libraries_cklaij/openmpi-1.4.5/bin/mpicxx<br>
Using Fortran linker: /opt/refresco/libraries_cklaij/openmpi-1.4.5/bin/mpif90<br>
Using libraries: -Wl,-rpath,/home/CKlaij/ReFRESCO/Dev/trunk/Libs/build/petsc-3.3-p1/linux_64bit_debug/lib -L/home/CKlaij/ReFRESCO/Dev/trunk/Libs/build/petsc-3.3-p1/linux_64bit_debug/lib -lpetsc -lX11 -lpthread -Wl,-rpath,/opt/refresco/libraries_cklaij/hypre-2.7.0b/lib -L/opt/refresco/libraries_cklaij/hypre-2.7.0b/lib -lHYPRE -Wl,-rpath,/opt/refresco/libraries_cklaij/ml-6.2/lib -L/opt/refresco/libraries_cklaij/ml-6.2/lib -lml -Wl,-rpath,/opt/intel/mkl -L/opt/intel/mkl -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -L/opt/refresco/libraries_cklaij/openmpi-1.4.5/lib -L/opt/intel/composer_xe_2011_sp1.9.293/compiler/lib/intel64 -L/opt/intel/composer_xe_2011_sp1.9.293/ipp/lib/intel64 -L/opt/intel/composer_xe_2011_sp1.9.293/mkl/lib/intel64 -L/opt/intel/composer_xe_2011_sp1.9.293/tbb/lib/intel64/cc4.1.0_libc2.4_kernel2.6.16.21 -L/usr/lib/gcc/x86_64-linux-gnu/4.4.3 -L/usr/lib/x86_64-linux-gnu -lmpi_f90 -lmpi_f77 -lifport -lifcore -lm -lm -lmpi_cxx -ldl -lmpi -lopen-rte -lopen-pal -lnsl -lutil -limf -lsvml -lipgo -ldecimal -lcilkrts -lstdc++ -lgcc_s -lirc -lpthread -lirc_s -ldl<br>
-----------------------------------------<br>
<br>
<br>
dr. ir. Christiaan Klaij<br>
CFD Researcher<br>
Research & Development<br>
E mailto:<a href="mailto:C.Klaij@marin.nl">C.Klaij@marin.nl</a><br>
T <a href="tel:%2B31%20317%2049%2033%2044" value="+31317493344">+31 317 49 33 44</a><br>
<br>
MARIN<br>
2, Haagsteeg, P.O. Box 28, 6700 AA Wageningen, The Netherlands<br>
T <a href="tel:%2B31%20317%2049%2039%2011" value="+31317493911">+31 317 49 39 11</a>, F <a href="tel:%2B31%20317%2049%2032%2045" value="+31317493245">+31 317 49 32 45</a>, I <a href="http://www.marin.nl" target="_blank">www.marin.nl</a><br>
<br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener<br>