May 2012 - petsc-dev - mailman3.cels.anl.gov

Re: [petsc-dev] API changes in MatIS
by Jed Brown 08 May '12

08 May '12

On Tue, May 8, 2012 at 7:32 AM, Stefano Zampini <stefano.zampini(a)gmail.com>wrote: > MatSetValues() is fine since only the global ordering is used. >> MatSetValuesLocal() would ideally continue to work with the old local space >> because the user has to redistribute their mesh to use the new ordering >> (maybe there should be a translation filter). >> >> > What do you mean for "translation filter"? > Translating a set of local indices for a different partition, but I misunderstood what you were suggesting. > > >> But how does this function actually work? Can the subdomain boundaries >> move, or are subdomains just moved to a different process? >> >> > Asking the first question: > The function works by first reordering local nodes with INTERNAL and > INTERFACE nodes (if the second argument LocalIS is not PETSC_NULL), in > order to be able to use memcpy when using local scatters of PCIS (e.g. > pcis->N_to_B). Then, it reorders the global numbering (if GlobalIS is not > PETSC_NULL) by exploiting data locality in the global vector and reduce the > amount of communications when scattering from global to local and > viceversa; in practice, a numbering of global nodes is constructed in such > a way that the greatest number of subdomain INTERNAL nodes will belong to > the local part of the global vector. Then, remaining INTERNAL nodes (notes > that if the parallel problem is well balanced, this will not occur) and > INTERFACE nodes are placed in the remaining global nodes, keeping them > sorted as they were in the previous global numbering. I think the latter > part can be further modified. When I'll add to matis.c, maybe you can take > a look at the code and give me your suggestions. > Okay, this is different from what I thought you were describing. Moving the global dofs is a tricky thing because the vectors we want to solve against and interpret results for come from a higher level. Changing the global ordering really produces a different matrix (MatMult is different), so it's not "optimization". If the linear operator stayed the same, but you just changed the format, perhaps we could call it "optimize". Are you aware of PCRedistribute? It could be made to perform the transformation you are thinking about for MATIS, solving the problem in the "optimized" ordering, and mapping the solution back to the original ordering. But why is this needed in the first place? Don't people who use MATIS also create a reasonable global ordering at setup time? It's bad for other parts of the algorithm to have a global numbering that doesn't line up well with the partition of elements. > > Regarding the other question: > With MATIS (and in principle every subassembled matrix) you are forced to > use the local matrices assembled by the user or an assembled group of them > (as you surely know since you are developing a PA). You cannot safely (in a > theoretical way) move subdomain boundaries since you don't know (at least > in principle) how to reproduce the matrix values on the local interface > since they are neumann boundaries for the local PDE. The main question is: > how can I split the assembled values on the nodes which previously were > internal to the subdomain? > This can't be done, at least without more information or making assumptions about the discretization. That's why I asked. > (I think this is the main question when thinking to a MatConvert > implementation which converts from a matrix obtained by DMDAGetMatrix to a > MATIS). > There is no reason why DMCreateMatrix() cannot return a MATIS. In any case, that's a better model than assembling in one format and trying to convert, especially since conversion really can't done in this direction (a MATIS can be converted to MATAIJ, not the other way). > > >> >> For a different, but related algorithm, I'm considering writing a MatPA >> (PA = partially assembled) which would choose some points to assemble >> internally. The reason for storing in partially assembled form is that I >> will want to solve a number of pinned Neumann problems to determine the >> coarse space and I'd rather not have to take unnecessary submatrices. Do >> you think this is a bad idea and I should just use MatIS? (I also want >> support for multiple subdomains per process and multiple processes per >> subdomain, though that could be added to MatIS.) >> > > Cool. Are you implementing FETI-DP? I'm planning to implement it in the > very next future for isogeometric analysis. Asking to your question, it > depends on how you need to solve these "pinned Neumann problems": you need > to solve unrelated local problems defined by the local matrices or solve a > parallel Neumann problem? In the first case, MatIS is your choice. In the > other case, I think you need to resort to a multigrid for a new MATIS > object defined by the subdomains of interest. > However, I will be glad to help you with PA if you want. I'm not that interested in FETI-DP since it's basically the same as BDDC and doesn't fix the exponential dependence on the number of levels. Also, I hate that it doesn't work well for irregular coefficients on subdomain boundaries. So I have a different algorithm in mind, more naturally viewed as a multigrid or multilevel Schwarz, but with aggressive coarsening enabled by the subdomains. One part of the method is compatible relaxation, for which pinned Neumann problems are used. The idea is to start with just enough pinning to make the subdomains non-singular, then enrich the coarse space by adaptively selecting vertices or edge/face moments that are not well controlled (as determined by compatible relaxation). Once the coarse injection is found, we construct the coarse basis functions by harmonic extension (again using pinned Neumann problems). With the multigrid hierarchy is constructed, the smoother would be based on overlapping subdomains or a multi-stage smoother involving both pinned Neumann problems and surface strips (or perhaps just point-block Jacobi). The problem with only using the pinned Neumann problems in the smoother seems to be that we can't remove the exponential dependence without allowing the smoother to move the injection (compatible relaxation is good for telling us whether a coarse space is good, but not actually the right thing for multigrid convergence).

1 0

Re: [petsc-dev] API changes in MatIS
by Stefano Zampini 08 May '12

08 May '12

> Okay, I just hate adding more places where such unnecessary constraints > are codified. > > Me too. I will take a look where to change the code to safely remove the requirement. Since we are speaking about MATIS. In my codes, I have a function (now it >> is in Fortran but I can translate it) >> >> MatISOptimize(Mat A,IS *LocalIS ,IS *GlobalIS) >> >> which changes the MATIS object, changing its local to global mapping to >> optimize either for local scatters and for global communications costs. It >> also changes the local matrix associated to the MATIS object. The function >> also returns the permutations used (if not PETSC_NULL). Can I add it to >> matis.c? Since it changes the underlying object, what are the requirements >> of PETSc? Should the user be able to insert the values with the original >> ordering? Any suggestions? > > > MatSetValues() is fine since only the global ordering is used. > MatSetValuesLocal() would ideally continue to work with the old local space > because the user has to redistribute their mesh to use the new ordering > (maybe there should be a translation filter). > > What do you mean for "translation filter"? > But how does this function actually work? Can the subdomain boundaries > move, or are subdomains just moved to a different process? > > Asking the first question: The function works by first reordering local nodes with INTERNAL and INTERFACE nodes (if the second argument LocalIS is not PETSC_NULL), in order to be able to use memcpy when using local scatters of PCIS (e.g. pcis->N_to_B). Then, it reorders the global numbering (if GlobalIS is not PETSC_NULL) by exploiting data locality in the global vector and reduce the amount of communications when scattering from global to local and viceversa; in practice, a numbering of global nodes is constructed in such a way that the greatest number of subdomain INTERNAL nodes will belong to the local part of the global vector. Then, remaining INTERNAL nodes (notes that if the parallel problem is well balanced, this will not occur) and INTERFACE nodes are placed in the remaining global nodes, keeping them sorted as they were in the previous global numbering. I think the latter part can be further modified. When I'll add to matis.c, maybe you can take a look at the code and give me your suggestions. Regarding the other question: With MATIS (and in principle every subassembled matrix) you are forced to use the local matrices assembled by the user or an assembled group of them (as you surely know since you are developing a PA). You cannot safely (in a theoretical way) move subdomain boundaries since you don't know (at least in principle) how to reproduce the matrix values on the local interface since they are neumann boundaries for the local PDE. The main question is: how can I split the assembled values on the nodes which previously were internal to the subdomain? (I think this is the main question when thinking to a MatConvert implementation which converts from a matrix obtained by DMDAGetMatrix to a MATIS). > > For a different, but related algorithm, I'm considering writing a MatPA > (PA = partially assembled) which would choose some points to assemble > internally. The reason for storing in partially assembled form is that I > will want to solve a number of pinned Neumann problems to determine the > coarse space and I'd rather not have to take unnecessary submatrices. Do > you think this is a bad idea and I should just use MatIS? (I also want > support for multiple subdomains per process and multiple processes per > subdomain, though that could be added to MatIS.) > Cool. Are you implementing FETI-DP? I'm planning to implement it in the very next future for isogeometric analysis. Asking to your question, it depends on how you need to solve these "pinned Neumann problems": you need to solve unrelated local problems defined by the local matrices or solve a parallel Neumann problem? In the first case, MatIS is your choice. In the other case, I think you need to resort to a multigrid for a new MATIS object defined by the subdomains of interest. However, I will be glad to help you with PA if you want. -- Stefano

1 0

Re: [petsc-dev] API changes in MatIS
by Jed Brown 08 May '12

08 May '12

On Tue, May 8, 2012 at 3:30 AM, Stefano Zampini <stefano.zampini(a)gmail.com>wrote: > Sorry, I didn't know that. A stupid question but, can you say how many > users are using MATIS objects? > I don't know, probably not a great number, but enough that it has to be documented. > > 2. What was wrong with the default MatGetVecs()? >> >> > There's a conceptual difference: default MatGetVecs uses mat->rmap->n and > mat->cmat->n as local sizes of the vectors since it assumes a distributed > Matrix. MATIS object is not distributed in a standard way (say, all > mat->rmap->n values sums to number of global rows), MatMult_IS is not > directly performed using the local part of the global vector, thus there is > not a real need (apart of Petsc checking the sizes) of the notion of > mat->rmap->n (or mat->cmap->n) as in the standard implementation. What you > really need is the global size of the vector; indeed the sizes of the local > vectors used to perform MatMult_IS doesn't sum up to the vector global size. > PetscLayout and the matrix size statements don't work as you describe here; if they do somewhere, it's broken. The matrix sizes are really statements about the distribution of the _vector_ that the matrix multiplies against and produces. With *AIJ, that happens to also be the row partition, and the off-diagonal part is managed in a compatible way, but that is not required. For dense matrices, it might make sense to store the matrix with an entirely different distribution (e.g. cyclic). An implementation issue like the sizes of the local matrices in MATIS should _not_ change the global matrix size (rmap and cmap). > > 3. Options prefixes should inherit the prefix from the parent and should >> end with "_" so that they aren't squashed together with the suboptions. >> > > Since MatSetOptionsPrefix(is->A,"is") was already there, I thought it was > right using the same prefix for local vectors. Change it as you prefer. > Ah, okay. It should be done like this in the future. http://petsc.cs.iit.edu/petsc/petsc-dev/rev/b70fc058783a > > >> 5. Why is the matrix being forced to be square? >> >> > MATIS were already forced to be square before I began using it. > Okay, I just hate adding more places where such unnecessary constraints are codified. > Pushed a fix to the problems you highlighted > > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/93e67397a066 > > Calling MatSetBlockSize before MatSetSizes in MatCreateIS solved all my > problems. > Great, thanks. > > Since we are speaking about MATIS. In my codes, I have a function (now it > is in Fortran but I can translate it) > > MatISOptimize(Mat A,IS *LocalIS ,IS *GlobalIS) > > which changes the MATIS object, changing its local to global mapping to > optimize either for local scatters and for global communications costs. It > also changes the local matrix associated to the MATIS object. The function > also returns the permutations used (if not PETSC_NULL). Can I add it to > matis.c? Since it changes the underlying object, what are the requirements > of PETSc? Should the user be able to insert the values with the original > ordering? Any suggestions? MatSetValues() is fine since only the global ordering is used. MatSetValuesLocal() would ideally continue to work with the old local space because the user has to redistribute their mesh to use the new ordering (maybe there should be a translation filter). But how does this function actually work? Can the subdomain boundaries move, or are subdomains just moved to a different process? For a different, but related algorithm, I'm considering writing a MatPA (PA = partially assembled) which would choose some points to assemble internally. The reason for storing in partially assembled form is that I will want to solve a number of pinned Neumann problems to determine the coarse space and I'd rather not have to take unnecessary submatrices. Do you think this is a bad idea and I should just use MatIS? (I also want support for multiple subdomains per process and multiple processes per subdomain, though that could be added to MatIS.)

1 0

Re: [petsc-dev] API changes in MatIS
by Stefano Zampini 08 May '12

08 May '12

All public API changes have to be documented at > src/docs/website/documentation/changes/dev.html. > > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/a942f19a37c9 > > Sorry, I didn't know that. A stupid question but, can you say how many users are using MATIS objects? 2. What was wrong with the default MatGetVecs()? > > There's a conceptual difference: default MatGetVecs uses mat->rmap->n and mat->cmat->n as local sizes of the vectors since it assumes a distributed Matrix. MATIS object is not distributed in a standard way (say, all mat->rmap->n values sums to number of global rows), MatMult_IS is not directly performed using the local part of the global vector, thus there is not a real need (apart of Petsc checking the sizes) of the notion of mat->rmap->n (or mat->cmap->n) as in the standard implementation. What you really need is the global size of the vector; indeed the sizes of the local vectors used to perform MatMult_IS doesn't sum up to the vector global size. 3. Options prefixes should inherit the prefix from the parent and should > end with "_" so that they aren't squashed together with the suboptions. > Since MatSetOptionsPrefix(is->A,"is") was already there, I thought it was right using the same prefix for local vectors. Change it as you prefer. > 5. Why is the matrix being forced to be square? > > MATIS were already forced to be square before I began using it. Pushed a fix to the problems you highlighted http://petsc.cs.iit.edu/petsc/petsc-dev/rev/93e67397a066 Calling MatSetBlockSize before MatSetSizes in MatCreateIS solved all my problems. Since we are speaking about MATIS. In my codes, I have a function (now it is in Fortran but I can translate it) MatISOptimize(Mat A,IS *LocalIS ,IS *GlobalIS) which changes the MATIS object, changing its local to global mapping to optimize either for local scatters and for global communications costs. It also changes the local matrix associated to the MATIS object. The function also returns the permutations used (if not PETSC_NULL). Can I add it to matis.c? Since it changes the underlying object, what are the requirements of PETSc? Should the user be able to insert the values with the original ordering? Any suggestions? -- Stefano

1 0

Re: [petsc-dev] Error during KSPDestroy
by Alexander Grayver 08 May '12

08 May '12

On 08.05.2012 08:59, Alexander Grayver wrote: > Barry, > > Not it works. I mean NOW, sorry :) > Thanks everybody! > > On 07.05.2012 22:36, Barry Smith wrote: >> Alexander >> >> Satish and I have determined the problem (took some valgrind and >> debugger work). We were not allocating enough "workspace" for one of >> the work arrays passed to zgesvd(). We have fixed it in petsc-dev. >> You should be able to do a hg pull -u then recompile the libraries >> with make cmake then relink and run your example. >> >> Thank you for your patience. >> >> Barry >> >> On May 7, 2012, at 9:10 AM, Alexander Grayver wrote: >> >>> On 07.05.2012 15:04, Barry Smith wrote: >>>> I am also running complex. >>>> >>>> Look in the file dlasq2.f (it will be in the externalpackages >>>> subdirectory of the PETSc directory. Look at line 215, this is >>>> where valgrind has a problem. In my copy >>>> >>>> END IF >>>> * >>>> * Check for negative data and compute sums of q's and e's. >>>> *<------ this is line 215 >>>> Z( 2*N ) = ZERO >>>> >>>> it is a comment, which is not good. Is lione 215 also a comment in >>>> your copy of dlasq2.f? >>> Barry, >>> >>> * >>> * Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...). >>> * >>> DO 30 K = 2*N, 2, -2 >>> Z( 2*K ) = ZERO >>> !<----------- LINE 215 >>> Z( 2*K-1 ) = Z( K ) >>> Z( 2*K-2 ) = ZERO >>> Z( 2*K-3 ) = Z( K-1 ) >>> 30 CONTINUE >>> >>> In valgrind log you can see that it complaints about following lines >>> as well: >>> >>> ==9009== Invalid write of size 8 >>> ==9009== at 0x10651D5: dlasq2_ (dlasq2.f:215) >>> ==9009== by 0x1064683: dlasq1_ (dlasq1.f:135) >>> ==9009== by 0x104EB3F: zbdsqr_ (zbdsqr.f:225) >>> ==9009== by 0x1023B74: zgesvd_ (zgesvd.f:2040) >>> ==9009== by 0xD38725: KSPComputeExtremeSingularValues_GMRES >>> (gmreig.c:46) >>> ==9009== by 0xCB3CC7: KSPComputeExtremeSingularValues (itfunc.c:47) >>> ==9009== by 0x406DF2: main (solveTest.c:47) >>> ==9009== Address 0x6ef5d88 is 8 bytes before a block of size 832 >>> alloc'd >>> ==9009== at 0x4C2786E: memalign (vg_replace_malloc.c:581) >>> ==9009== by 0x47E3CB: PetscMallocAlign (mal.c:30) >>> ==9009== by 0xD2E286: KSPSetUp_GMRES (gmres.c:73) >>> ==9009== by 0xCB5464: KSPSetUp (itfunc.c:239) >>> ==9009== by 0xCB6E56: KSPSolve (itfunc.c:402) >>> ==9009== by 0x406DDB: main (solveTest.c:46) >>> ==9009== >>> ==9009== Invalid write of size 8 >>> ==9009== at 0x1065204: dlasq2_ (dlasq2.f:216) >>> .... >>> ==9009== >>> ==9009== Invalid write of size 8 >>> ==9009== at 0x1065223: dlasq2_ (dlasq2.f:217) >>> .... >>> ==9009== >>> ==9009== Invalid write of size 8 >>> ==9009== at 0x1065255: dlasq2_ (dlasq2.f:218) >>> .... >>> >>> All further output is also related to the Z array. >>> Hard to believe this is a LAPACK problem... I tried 3 >>> implementations over 2 machines. >>> I have bad feeling it's my stupid mistake somewhere... :) >>> >>> Just in case, I run ubuntu 11.1 and PETSc is configured like this >>> with default gcc compiler: >>> ./configure --with-petsc-arch=mpich-gcc-complex-debug-c >>> --download-f-blas-lapack --with-precision=double >>> --with-scalar-type=complex --download-mpich >>> >>>> There are two possible causes I can think of for your problem >>>> >>>> 1) PETSc does not allocate enough work space for zgesvd() or >>>> 2) the BLAS/LAPACK routines have a bug where they sometimes access >>>> out of their work space. >>>> >>>> >>>> Satish, >>>> >>>> Can you try the same build options on a Linux machine as >>>> close to Alexander as we have and see if you can reproduce this? >>>> >>>> >>>> Barry >>>> >>>> >>>> >>>> On May 7, 2012, at 2:16 AM, Alexander Grayver wrote: >>>> >>>>> On 06.05.2012 22:24, Barry Smith wrote: >>>>>> Alexander, >>>>>> >>>>>> I cannot reproduce this on my mac with 3 different >>>>>> blas/lapack. >>>>> Barry, >>>>> >>>>> I'm surprised. I ran it on my home PC with ubuntu and PETSc >>>>> configured from scratch as following: >>>>> --download-mpich --with-fortran-interfaces=1 --download-scalapack >>>>> --download-blacs --with-scalar-type=complex --download-blas-lapack >>>>> --with-precision=double >>>>> >>>>> And it's still there. >>>>> Please note that all my numbers are complex. >>>>> >>>>>> Could you please run the case below but with >>>>>> --download-f-blas-lapack (you forgot the -f last time)? Send us >>>>>> the valgrind results. This will tell use the exact line number in >>>>>> dlasq3() that is triggering the bad read. >>>>> I did: >>>>> ./configure --with-petsc-arch=openmpi-intel-complex-debug-c >>>>> --download-scalapack --download-blacs --download-f-blas-lapack >>>>> --with-precision=double --with-scalar-type=complex >>>>> >>>>> And then valgrind program. The first message from log: >>>>> >>>>> ==27656== Invalid write of size 8 >>>>> ==27656== at 0x15A8E9E: dlasq2_ (dlasq2.f:215) >>>>> ==27656== by 0x15A83A4: dlasq1_ (dlasq1.f:135) >>>>> ==27656== by 0x158ACEC: zbdsqr_ (zbdsqr.f:225) >>>>> ==27656== by 0x154EC27: zgesvd_ (zgesvd.f:2038) >>>>> ==27656== by 0x695DD3: KSPComputeExtremeSingularValues_GMRES >>>>> (gmreig.c:46) >>>>> ==27656== by 0x69DD76: KSPComputeExtremeSingularValues >>>>> (itfunc.c:47) >>>>> ==27656== by 0x44E98C: main (solveTest.c:62) >>>>> ==27656== Address 0xfad2d98 is 8 bytes before a block of size 832 >>>>> alloc'd >>>>> ==27656== at 0x4C25D66: memalign (vg_replace_malloc.c:694) >>>>> ==27656== by 0x4B642B: PetscMallocAlign (mal.c:30) >>>>> ==27656== by 0x687775: KSPSetUp_GMRES (gmres.c:73) >>>>> ==27656== by 0x69FE4A: KSPSetUp (itfunc.c:239) >>>>> ==27656== by 0x6A2058: KSPSolve (itfunc.c:402) >>>>> ==27656== by 0x44E969: main (solveTest.c:61) >>>>> >>>>> Please find full log attached. >>>>> >>>>>> Barry >>>>>> >>>>>> >>>>>> On May 6, 2012, at 9:16 AM, Alexander Grayver wrote: >>>>>> >>>>>>> On 06.05.2012 15:34, Matthew Knepley wrote: >>>>>>>> On Sun, May 6, 2012 at 9:24 AM, Alexander >>>>>>>> Grayver<agrayver(a)gfz-potsdam.de> wrote: >>>>>>>> Hm, valgrind gives a lot of output like that (see full log in >>>>>>>> previous message): >>>>>>>> >>>>>>>> Can you run this with --download-f-blas-lapack? This sounds >>>>>>>> much more like an MKL bug. >>>>>>> I did: >>>>>>> --download-scalapack --download-blacs --download-blas-lapack >>>>>>> --with-precision=double --with-scalar-type=complex >>>>>>> >>>>>>> The error is still there. I checked "ldd solveTest", mkl is not >>>>>>> used for sure. This is not an MKL problem I guess: >>>>>>> >>>>>>> ==13600== Invalid read of size 8 >>>>>>> ==13600== at 0x58636AF: dlasq3_ (in >>>>>>> /usr/local/lib/liblapack.so.3.2.2) >>>>>>> ==13600== by 0x5862C84: dlasq2_ (in >>>>>>> /usr/local/lib/liblapack.so.3.2.2) >>>>>>> ==13600== by 0x5861F2C: dlasq1_ (in >>>>>>> /usr/local/lib/liblapack.so.3.2.2) >>>>>>> ==13600== by 0x571A479: zbdsqr_ (in >>>>>>> /usr/local/lib/liblapack.so.3.2.2) >>>>>>> ==13600== by 0x57466A7: zgesvd_ (in >>>>>>> /usr/local/lib/liblapack.so.3.2.2) >>>>>>> ==13600== by 0x694687: KSPComputeExtremeSingularValues_GMRES >>>>>>> (gmreig.c:46) >>>>>>> ==13600== by 0x69C62A: KSPComputeExtremeSingularValues >>>>>>> (itfunc.c:47) >>>>>>> ==13600== by 0x44E02C: main (solveTest.c:62) >>>>>>> ==13600== Address 0x10826b90 is 16 bytes before a block of size >>>>>>> 832 alloc'd >>>>>>> ==13600== at 0x4C25D66: memalign (vg_replace_malloc.c:694) >>>>>>> ==13600== by 0x4B5ACB: PetscMallocAlign (mal.c:30) >>>>>>> ==13600== by 0x686181: KSPSetUp_GMRES (gmres.c:73) >>>>>>> ==13600== by 0x69E6FE: KSPSetUp (itfunc.c:239) >>>>>>> ==13600== by 0x6A090C: KSPSolve (itfunc.c:402) >>>>>>> ==13600== by 0x44E009: main (solveTest.c:61) >>>>>>> >>>>>>> The weird thing is that the it gives correct result, so zgesvd >>>>>>> works fine. >>>>>>> >>>>>>> And also running this program with 10 iterations in valgrind >>>>>>> doesn't produce error. The low above is with 100 iterations. >>>>>>> Without valgrind the error is always there. >>>>>>> >>>>>>> -- >>>>>>> Regards, >>>>>>> Alexander >>>>>>> >>>>> -- >>>>> Regards, >>>>> Alexander >>>>> >>>>> <valgrind.zip> >>> >>> -- >>> Regards, >>> Alexander >>> > > -- Regards, Alexander

1 0

Re: [petsc-dev] Error during KSPDestroy
by Alexander Grayver 08 May '12

08 May '12

Barry, Not it works. Thanks everybody! On 07.05.2012 22:36, Barry Smith wrote: > Alexander > > Satish and I have determined the problem (took some valgrind and debugger work). We were not allocating enough "workspace" for one of the work arrays passed to zgesvd(). We have fixed it in petsc-dev. You should be able to do a hg pull -u then recompile the libraries with make cmake then relink and run your example. > > Thank you for your patience. > > Barry > > On May 7, 2012, at 9:10 AM, Alexander Grayver wrote: > >> On 07.05.2012 15:04, Barry Smith wrote: >>> I am also running complex. >>> >>> Look in the file dlasq2.f (it will be in the externalpackages subdirectory of the PETSc directory. Look at line 215, this is where valgrind has a problem. In my copy >>> >>> END IF >>> * >>> * Check for negative data and compute sums of q's and e's. >>> *<------ this is line 215 >>> Z( 2*N ) = ZERO >>> >>> it is a comment, which is not good. Is lione 215 also a comment in your copy of dlasq2.f? >> Barry, >> >> * >> * Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...). >> * >> DO 30 K = 2*N, 2, -2 >> Z( 2*K ) = ZERO !<----------- LINE 215 >> Z( 2*K-1 ) = Z( K ) >> Z( 2*K-2 ) = ZERO >> Z( 2*K-3 ) = Z( K-1 ) >> 30 CONTINUE >> >> In valgrind log you can see that it complaints about following lines as well: >> >> ==9009== Invalid write of size 8 >> ==9009== at 0x10651D5: dlasq2_ (dlasq2.f:215) >> ==9009== by 0x1064683: dlasq1_ (dlasq1.f:135) >> ==9009== by 0x104EB3F: zbdsqr_ (zbdsqr.f:225) >> ==9009== by 0x1023B74: zgesvd_ (zgesvd.f:2040) >> ==9009== by 0xD38725: KSPComputeExtremeSingularValues_GMRES (gmreig.c:46) >> ==9009== by 0xCB3CC7: KSPComputeExtremeSingularValues (itfunc.c:47) >> ==9009== by 0x406DF2: main (solveTest.c:47) >> ==9009== Address 0x6ef5d88 is 8 bytes before a block of size 832 alloc'd >> ==9009== at 0x4C2786E: memalign (vg_replace_malloc.c:581) >> ==9009== by 0x47E3CB: PetscMallocAlign (mal.c:30) >> ==9009== by 0xD2E286: KSPSetUp_GMRES (gmres.c:73) >> ==9009== by 0xCB5464: KSPSetUp (itfunc.c:239) >> ==9009== by 0xCB6E56: KSPSolve (itfunc.c:402) >> ==9009== by 0x406DDB: main (solveTest.c:46) >> ==9009== >> ==9009== Invalid write of size 8 >> ==9009== at 0x1065204: dlasq2_ (dlasq2.f:216) >> .... >> ==9009== >> ==9009== Invalid write of size 8 >> ==9009== at 0x1065223: dlasq2_ (dlasq2.f:217) >> .... >> ==9009== >> ==9009== Invalid write of size 8 >> ==9009== at 0x1065255: dlasq2_ (dlasq2.f:218) >> .... >> >> All further output is also related to the Z array. >> Hard to believe this is a LAPACK problem... I tried 3 implementations over 2 machines. >> I have bad feeling it's my stupid mistake somewhere... :) >> >> Just in case, I run ubuntu 11.1 and PETSc is configured like this with default gcc compiler: >> ./configure --with-petsc-arch=mpich-gcc-complex-debug-c --download-f-blas-lapack --with-precision=double --with-scalar-type=complex --download-mpich >> >>> There are two possible causes I can think of for your problem >>> >>> 1) PETSc does not allocate enough work space for zgesvd() or >>> 2) the BLAS/LAPACK routines have a bug where they sometimes access out of their work space. >>> >>> >>> Satish, >>> >>> Can you try the same build options on a Linux machine as close to Alexander as we have and see if you can reproduce this? >>> >>> >>> Barry >>> >>> >>> >>> On May 7, 2012, at 2:16 AM, Alexander Grayver wrote: >>> >>>> On 06.05.2012 22:24, Barry Smith wrote: >>>>> Alexander, >>>>> >>>>> I cannot reproduce this on my mac with 3 different blas/lapack. >>>> Barry, >>>> >>>> I'm surprised. I ran it on my home PC with ubuntu and PETSc configured from scratch as following: >>>> --download-mpich --with-fortran-interfaces=1 --download-scalapack --download-blacs --with-scalar-type=complex --download-blas-lapack --with-precision=double >>>> >>>> And it's still there. >>>> Please note that all my numbers are complex. >>>> >>>>> Could you please run the case below but with --download-f-blas-lapack (you forgot the -f last time)? Send us the valgrind results. This will tell use the exact line number in dlasq3() that is triggering the bad read. >>>> I did: >>>> ./configure --with-petsc-arch=openmpi-intel-complex-debug-c --download-scalapack --download-blacs --download-f-blas-lapack --with-precision=double --with-scalar-type=complex >>>> >>>> And then valgrind program. The first message from log: >>>> >>>> ==27656== Invalid write of size 8 >>>> ==27656== at 0x15A8E9E: dlasq2_ (dlasq2.f:215) >>>> ==27656== by 0x15A83A4: dlasq1_ (dlasq1.f:135) >>>> ==27656== by 0x158ACEC: zbdsqr_ (zbdsqr.f:225) >>>> ==27656== by 0x154EC27: zgesvd_ (zgesvd.f:2038) >>>> ==27656== by 0x695DD3: KSPComputeExtremeSingularValues_GMRES (gmreig.c:46) >>>> ==27656== by 0x69DD76: KSPComputeExtremeSingularValues (itfunc.c:47) >>>> ==27656== by 0x44E98C: main (solveTest.c:62) >>>> ==27656== Address 0xfad2d98 is 8 bytes before a block of size 832 alloc'd >>>> ==27656== at 0x4C25D66: memalign (vg_replace_malloc.c:694) >>>> ==27656== by 0x4B642B: PetscMallocAlign (mal.c:30) >>>> ==27656== by 0x687775: KSPSetUp_GMRES (gmres.c:73) >>>> ==27656== by 0x69FE4A: KSPSetUp (itfunc.c:239) >>>> ==27656== by 0x6A2058: KSPSolve (itfunc.c:402) >>>> ==27656== by 0x44E969: main (solveTest.c:61) >>>> >>>> Please find full log attached. >>>> >>>>> Barry >>>>> >>>>> >>>>> On May 6, 2012, at 9:16 AM, Alexander Grayver wrote: >>>>> >>>>>> On 06.05.2012 15:34, Matthew Knepley wrote: >>>>>>> On Sun, May 6, 2012 at 9:24 AM, Alexander Grayver<agrayver(a)gfz-potsdam.de> wrote: >>>>>>> Hm, valgrind gives a lot of output like that (see full log in previous message): >>>>>>> >>>>>>> Can you run this with --download-f-blas-lapack? This sounds much more like an MKL bug. >>>>>> I did: >>>>>> --download-scalapack --download-blacs --download-blas-lapack --with-precision=double --with-scalar-type=complex >>>>>> >>>>>> The error is still there. I checked "ldd solveTest", mkl is not used for sure. This is not an MKL problem I guess: >>>>>> >>>>>> ==13600== Invalid read of size 8 >>>>>> ==13600== at 0x58636AF: dlasq3_ (in /usr/local/lib/liblapack.so.3.2.2) >>>>>> ==13600== by 0x5862C84: dlasq2_ (in /usr/local/lib/liblapack.so.3.2.2) >>>>>> ==13600== by 0x5861F2C: dlasq1_ (in /usr/local/lib/liblapack.so.3.2.2) >>>>>> ==13600== by 0x571A479: zbdsqr_ (in /usr/local/lib/liblapack.so.3.2.2) >>>>>> ==13600== by 0x57466A7: zgesvd_ (in /usr/local/lib/liblapack.so.3.2.2) >>>>>> ==13600== by 0x694687: KSPComputeExtremeSingularValues_GMRES (gmreig.c:46) >>>>>> ==13600== by 0x69C62A: KSPComputeExtremeSingularValues (itfunc.c:47) >>>>>> ==13600== by 0x44E02C: main (solveTest.c:62) >>>>>> ==13600== Address 0x10826b90 is 16 bytes before a block of size 832 alloc'd >>>>>> ==13600== at 0x4C25D66: memalign (vg_replace_malloc.c:694) >>>>>> ==13600== by 0x4B5ACB: PetscMallocAlign (mal.c:30) >>>>>> ==13600== by 0x686181: KSPSetUp_GMRES (gmres.c:73) >>>>>> ==13600== by 0x69E6FE: KSPSetUp (itfunc.c:239) >>>>>> ==13600== by 0x6A090C: KSPSolve (itfunc.c:402) >>>>>> ==13600== by 0x44E009: main (solveTest.c:61) >>>>>> >>>>>> The weird thing is that the it gives correct result, so zgesvd works fine. >>>>>> >>>>>> And also running this program with 10 iterations in valgrind doesn't produce error. The low above is with 100 iterations. >>>>>> Without valgrind the error is always there. >>>>>> >>>>>> -- >>>>>> Regards, >>>>>> Alexander >>>>>> >>>> -- >>>> Regards, >>>> Alexander >>>> >>>> <valgrind.zip> >> >> -- >> Regards, >> Alexander >> -- Regards, Alexander

1 0

API changes in MatIS
by Jed Brown 07 May '12

07 May '12

All public API changes have to be documented at src/docs/website/documentation/changes/dev.html. http://petsc.cs.iit.edu/petsc/petsc-dev/rev/a942f19a37c9 Some questions and comments on this patch: 1. Why don't you use the API (e.g. PetscLayoutSetBlockSize()) when setting the block size? 2. What was wrong with the default MatGetVecs()? 3. Options prefixes should inherit the prefix from the parent and should end with "_" so that they aren't squashed together with the suboptions. 4. MatCreateIS() leaks the new Vec "global". 5. Why is the matrix being forced to be square? 6. Why is the size of the rmap/cmap being changed to be possibly different from m,n,M,N passed in? By not using the PetscLayout API, it looks like we can end up with a corrupt PetscLayout this way.

1 0

Re: [petsc-dev] Error during KSPDestroy
by Barry Smith 07 May '12

07 May '12

Alexander Satish and I have determined the problem (took some valgrind and debugger work). We were not allocating enough "workspace" for one of the work arrays passed to zgesvd(). We have fixed it in petsc-dev. You should be able to do a hg pull -u then recompile the libraries with make cmake then relink and run your example. Thank you for your patience. Barry On May 7, 2012, at 9:10 AM, Alexander Grayver wrote: > On 07.05.2012 15:04, Barry Smith wrote: >> I am also running complex. >> >> Look in the file dlasq2.f (it will be in the externalpackages subdirectory of the PETSc directory. Look at line 215, this is where valgrind has a problem. In my copy >> >> END IF >> * >> * Check for negative data and compute sums of q's and e's. >> *<------ this is line 215 >> Z( 2*N ) = ZERO >> >> it is a comment, which is not good. Is lione 215 also a comment in your copy of dlasq2.f? > Barry, > > * > * Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...). > * > DO 30 K = 2*N, 2, -2 > Z( 2*K ) = ZERO ! <----------- LINE 215 > Z( 2*K-1 ) = Z( K ) > Z( 2*K-2 ) = ZERO > Z( 2*K-3 ) = Z( K-1 ) > 30 CONTINUE > > In valgrind log you can see that it complaints about following lines as well: > > ==9009== Invalid write of size 8 > ==9009== at 0x10651D5: dlasq2_ (dlasq2.f:215) > ==9009== by 0x1064683: dlasq1_ (dlasq1.f:135) > ==9009== by 0x104EB3F: zbdsqr_ (zbdsqr.f:225) > ==9009== by 0x1023B74: zgesvd_ (zgesvd.f:2040) > ==9009== by 0xD38725: KSPComputeExtremeSingularValues_GMRES (gmreig.c:46) > ==9009== by 0xCB3CC7: KSPComputeExtremeSingularValues (itfunc.c:47) > ==9009== by 0x406DF2: main (solveTest.c:47) > ==9009== Address 0x6ef5d88 is 8 bytes before a block of size 832 alloc'd > ==9009== at 0x4C2786E: memalign (vg_replace_malloc.c:581) > ==9009== by 0x47E3CB: PetscMallocAlign (mal.c:30) > ==9009== by 0xD2E286: KSPSetUp_GMRES (gmres.c:73) > ==9009== by 0xCB5464: KSPSetUp (itfunc.c:239) > ==9009== by 0xCB6E56: KSPSolve (itfunc.c:402) > ==9009== by 0x406DDB: main (solveTest.c:46) > ==9009== > ==9009== Invalid write of size 8 > ==9009== at 0x1065204: dlasq2_ (dlasq2.f:216) > .... > ==9009== > ==9009== Invalid write of size 8 > ==9009== at 0x1065223: dlasq2_ (dlasq2.f:217) > .... > ==9009== > ==9009== Invalid write of size 8 > ==9009== at 0x1065255: dlasq2_ (dlasq2.f:218) > .... > > All further output is also related to the Z array. > Hard to believe this is a LAPACK problem... I tried 3 implementations over 2 machines. > I have bad feeling it's my stupid mistake somewhere... :) > > Just in case, I run ubuntu 11.1 and PETSc is configured like this with default gcc compiler: > ./configure --with-petsc-arch=mpich-gcc-complex-debug-c --download-f-blas-lapack --with-precision=double --with-scalar-type=complex --download-mpich > >> There are two possible causes I can think of for your problem >> >> 1) PETSc does not allocate enough work space for zgesvd() or >> 2) the BLAS/LAPACK routines have a bug where they sometimes access out of their work space. >> >> >> Satish, >> >> Can you try the same build options on a Linux machine as close to Alexander as we have and see if you can reproduce this? >> >> >> Barry >> >> >> >> On May 7, 2012, at 2:16 AM, Alexander Grayver wrote: >> >>> On 06.05.2012 22:24, Barry Smith wrote: >>>> Alexander, >>>> >>>> I cannot reproduce this on my mac with 3 different blas/lapack. >>> Barry, >>> >>> I'm surprised. I ran it on my home PC with ubuntu and PETSc configured from scratch as following: >>> --download-mpich --with-fortran-interfaces=1 --download-scalapack --download-blacs --with-scalar-type=complex --download-blas-lapack --with-precision=double >>> >>> And it's still there. >>> Please note that all my numbers are complex. >>> >>>> Could you please run the case below but with --download-f-blas-lapack (you forgot the -f last time)? Send us the valgrind results. This will tell use the exact line number in dlasq3() that is triggering the bad read. >>> I did: >>> ./configure --with-petsc-arch=openmpi-intel-complex-debug-c --download-scalapack --download-blacs --download-f-blas-lapack --with-precision=double --with-scalar-type=complex >>> >>> And then valgrind program. The first message from log: >>> >>> ==27656== Invalid write of size 8 >>> ==27656== at 0x15A8E9E: dlasq2_ (dlasq2.f:215) >>> ==27656== by 0x15A83A4: dlasq1_ (dlasq1.f:135) >>> ==27656== by 0x158ACEC: zbdsqr_ (zbdsqr.f:225) >>> ==27656== by 0x154EC27: zgesvd_ (zgesvd.f:2038) >>> ==27656== by 0x695DD3: KSPComputeExtremeSingularValues_GMRES (gmreig.c:46) >>> ==27656== by 0x69DD76: KSPComputeExtremeSingularValues (itfunc.c:47) >>> ==27656== by 0x44E98C: main (solveTest.c:62) >>> ==27656== Address 0xfad2d98 is 8 bytes before a block of size 832 alloc'd >>> ==27656== at 0x4C25D66: memalign (vg_replace_malloc.c:694) >>> ==27656== by 0x4B642B: PetscMallocAlign (mal.c:30) >>> ==27656== by 0x687775: KSPSetUp_GMRES (gmres.c:73) >>> ==27656== by 0x69FE4A: KSPSetUp (itfunc.c:239) >>> ==27656== by 0x6A2058: KSPSolve (itfunc.c:402) >>> ==27656== by 0x44E969: main (solveTest.c:61) >>> >>> Please find full log attached. >>> >>>> Barry >>>> >>>> >>>> On May 6, 2012, at 9:16 AM, Alexander Grayver wrote: >>>> >>>>> On 06.05.2012 15:34, Matthew Knepley wrote: >>>>>> On Sun, May 6, 2012 at 9:24 AM, Alexander Grayver<agrayver(a)gfz-potsdam.de> wrote: >>>>>> Hm, valgrind gives a lot of output like that (see full log in previous message): >>>>>> >>>>>> Can you run this with --download-f-blas-lapack? This sounds much more like an MKL bug. >>>>> I did: >>>>> --download-scalapack --download-blacs --download-blas-lapack --with-precision=double --with-scalar-type=complex >>>>> >>>>> The error is still there. I checked "ldd solveTest", mkl is not used for sure. This is not an MKL problem I guess: >>>>> >>>>> ==13600== Invalid read of size 8 >>>>> ==13600== at 0x58636AF: dlasq3_ (in /usr/local/lib/liblapack.so.3.2.2) >>>>> ==13600== by 0x5862C84: dlasq2_ (in /usr/local/lib/liblapack.so.3.2.2) >>>>> ==13600== by 0x5861F2C: dlasq1_ (in /usr/local/lib/liblapack.so.3.2.2) >>>>> ==13600== by 0x571A479: zbdsqr_ (in /usr/local/lib/liblapack.so.3.2.2) >>>>> ==13600== by 0x57466A7: zgesvd_ (in /usr/local/lib/liblapack.so.3.2.2) >>>>> ==13600== by 0x694687: KSPComputeExtremeSingularValues_GMRES (gmreig.c:46) >>>>> ==13600== by 0x69C62A: KSPComputeExtremeSingularValues (itfunc.c:47) >>>>> ==13600== by 0x44E02C: main (solveTest.c:62) >>>>> ==13600== Address 0x10826b90 is 16 bytes before a block of size 832 alloc'd >>>>> ==13600== at 0x4C25D66: memalign (vg_replace_malloc.c:694) >>>>> ==13600== by 0x4B5ACB: PetscMallocAlign (mal.c:30) >>>>> ==13600== by 0x686181: KSPSetUp_GMRES (gmres.c:73) >>>>> ==13600== by 0x69E6FE: KSPSetUp (itfunc.c:239) >>>>> ==13600== by 0x6A090C: KSPSolve (itfunc.c:402) >>>>> ==13600== by 0x44E009: main (solveTest.c:61) >>>>> >>>>> The weird thing is that the it gives correct result, so zgesvd works fine. >>>>> >>>>> And also running this program with 10 iterations in valgrind doesn't produce error. The low above is with 100 iterations. >>>>> Without valgrind the error is always there. >>>>> >>>>> -- >>>>> Regards, >>>>> Alexander >>>>> >>> >>> -- >>> Regards, >>> Alexander >>> >>> <valgrind.zip> > > > -- > Regards, > Alexander >

1 0

Re: [petsc-dev] Error during KSPDestroy
by Barry Smith 07 May '12

07 May '12

On May 7, 2012, at 11:18 AM, Matthew Knepley wrote: > On Mon, May 7, 2012 at 11:53 AM, Satish Balay <balay(a)mcs.anl.gov> wrote: > On Mon, 7 May 2012, Barry Smith wrote: > > > > > I am also running complex. > > > > Look in the file dlasq2.f (it will be in the externalpackages subdirectory of the PETSc directory. Look at line 215, this is where valgrind has a problem. In my copy > > > > END IF > > * > > * Check for negative data and compute sums of q's and e's. > > * <------ this is line 215 > > Z( 2*N ) = ZERO > > > > it is a comment, which is not good. Is lione 215 also a comment in your copy of dlasq2.f? > > > > There are two possible causes I can think of for your problem > > > > 1) PETSc does not allocate enough work space for zgesvd() or > > 2) the BLAS/LAPACK routines have a bug where they sometimes access out of their work space. > > > > > > Satish, > > > > Can you try the same build options on a Linux machine as close to Alexander as we have and see if you can reproduce this? > > I can reproduce on my linux laptop. [also on petsc.mcs] > > Are we forgetting the factor 2 on the workspace? No because we are pure complex so all our mallocs are with respect to complex. I looked through the docs for the lapack routine and it looks like we are allocating enough. But. Satish, Can you bring your laptop up and we'll see if we can track down the problem? Thanks Barry > > Matt > > > Satish > -- > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > -- Norbert Wiener

1 0

Re: [petsc-dev] Error during KSPDestroy
by Matthew Knepley 07 May '12

07 May '12

On Mon, May 7, 2012 at 11:53 AM, Satish Balay <balay(a)mcs.anl.gov> wrote: > On Mon, 7 May 2012, Barry Smith wrote: > > > > > I am also running complex. > > > > Look in the file dlasq2.f (it will be in the externalpackages > subdirectory of the PETSc directory. Look at line 215, this is where > valgrind has a problem. In my copy > > > > END IF > > * > > * Check for negative data and compute sums of q's and e's. > > * > <------ this is line 215 > > Z( 2*N ) = ZERO > > > > it is a comment, which is not good. Is lione 215 also a comment in your > copy of dlasq2.f? > > > > There are two possible causes I can think of for your problem > > > > 1) PETSc does not allocate enough work space for zgesvd() or > > 2) the BLAS/LAPACK routines have a bug where they sometimes access out > of their work space. > > > > > > Satish, > > > > Can you try the same build options on a Linux machine as close to > Alexander as we have and see if you can reproduce this? > > I can reproduce on my linux laptop. [also on petsc.mcs] Are we forgetting the factor 2 on the workspace? Matt > > Satish -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener

1 0