Re: [petsc-users] Strange efficiency in PETSc-dev using OpenMP
I think this is definitely an issue with setting the affinities for threads, i.e., the assignment of threads to cores. Ideally each thread should be assigned to a distinct core but in your case all the 4 threads are getting pinned to the same core resulting in such a massive slowdown. Unfortunately, the thread affinities for OpenMP are set through environment variables. For Intel's OpenMP one needs to define the thread affinities through the environment variable KMP_AFFINITY. See this document here http://software.intel.com/sites/products/documentation/studio/composer/en-us/2011Update/compiler_c/optaps/common/optaps_openmp_thread_affinity.htm. Try setting the affinities via KMP_AFFINITY and let us know if it works. Shri On Sep 21, 2013, at 11:06 PM, Danyang Su wrote: Hi Shri, Thanks for your info. It can work with the option -threadcomm_type openmp. But another problem arises, as described as follows. The sparse matrix is 53760*53760 with 1067392 non-zero entries. If the codes is compiled using PETSc-3.4.2, it works fine, the equations can be solved quickly and I can see the speedup. But if the code is compiled using PETSc-dev with OpenMP option, it takes a long time in solving the equations and I cannot see any speedup when more processors are used. For PETSc-3.4.2, run by mpiexec -n 4 ksp_inhm_d -log_summary log_mpi4_petsc3.4.2.log, the iteration and runtime are: Iterations 6 time_assembly 0.4137E-01 time_ksp 0.9296E-01 For PETSc-dev, run by mpiexec -n 1 ksp_inhm_d -threadcomm_type openmp -threadcomm_nthreads 4 -log_summary log_openmp_petsc_dev.log, the iteration and runtime are: Iterations 6 time_assembly 0.3595E+03 time_ksp 0.2907E+00 Most of the time 'time_assembly 0.3595E+03' is spent on the following codes do i = istart, iend - 1 ii = ia_in(i+1) jj = ia_in(i+2) call MatSetValues(a, ione, i, jj-ii, ja_in(ii:jj-1)-1, a_in(ii:jj-1), Insert_Values, ierr) end do The log files for both PETSc-3.4.2 and PETSc-dev are attached. Is there anything wrong with my codes or with running option? The above codes works fine when using MPICH. Thanks and regards, Danyang On 21/09/2013 2:09 PM, Shri wrote: There are three thread communicator types in PETSc. The default is no thread which is basically a non-threaded version. The other two types are openmp and pthread. If you want to use OpenMP then use the option -threadcomm_type openmp. Shri On Sep 21, 2013, at 3:46 PM, Danyang Su danyang...@gmail.com wrote: Hi Barry, Thanks for the quick reply. After changing #if defined(PETSC_HAVE_PTHREADCLASSES) || defined (PETSC_HAVE_OPENMP) to #if defined(PETSC_HAVE_PTHREADCLASSES) and comment out #elif defined(PETSC_HAVE_OPENMP) PETSC_EXTERN PetscStack *petscstack; It can be compiled and validated with make test. But I still have questions on running the examples. After rebuild the codes (e.g., ksp_ex2f.f), I can run it with mpiexec -n 1 ksp_ex2f, or mpiexec -n 4 ksp_ex2f, or mpiexec -n 1 ksp_ex2f -threadcomm_nthreads 1, but if I run it with mpiexec -n 1 ksp_ex2f -threadcomm_nthreads 4, there will be a lot of error information (attached). The codes is not modified and there is no OpenMP routines in it. For the current development in my project, I want to keep the OpenMP codes in calculating matrix values, but want to solve it with PETSc (OpenMP). Is it possible? Thanks and regards, Danyang On 21/09/2013 7:26 AM, Barry Smith wrote: Danyang, I don't think the || defined (PETSC_HAVE_OPENMP) belongs in the code below. /* Linux functions CPU_SET and others don't work if sched.h is not included before including pthread.h. Also, these functions are active only if either _GNU_SOURCE or __USE_GNU is not set (see /usr/include/sched.h and /usr/include/features.h), hence set these first. */ #if defined(PETSC_HAVE_PTHREADCLASSES) || defined (PETSC_HAVE_OPENMP) Edit include/petscerror.h and locate these lines and remove that part and then rerun make all. Let us know if it works or not. Barry i.e. replace #if defined(PETSC_HAVE_PTHREADCLASSES) || defined (PETSC_HAVE_OPENMP) with #if defined(PETSC_HAVE_PTHREADCLASSES) On Sep 21, 2013, at 6:53 AM, Matthew Knepley petsc-ma...@mcs.anl.gov wrote: On Sat, Sep 21, 2013 at 12:18 AM, Danyang Su danyang...@gmail.com wrote: Hi All, I got error information in compiling petsc-dev with openmp in cygwin. Before, I have successfully compiled petsc-3.4.2 and it works fine. The log files have been attached. The OpenMP configure test is wrong. It clearly fails to find pthread.h, but the test passes. Then in petscerror.h we guard pthread.h using PETSC_HAVE_OPENMP. Can someone who knows OpenMP fix this? Matt Thanks, Danyang -- What most experimenters take
Re: [petsc-users] Performance of PETSc TS solver
Is it possible for you to share the IJacobian code? That would help us to understand the issue better and faster. Thanks, Shri On Aug 16, 2013, at 4:07 PM, Barry Smith wrote: On Aug 16, 2013, at 4:01 PM, Jin, Shuangshuang shuangshuang@pnnl.gov wrote: Hello, Jed, in my IJacobian subroutine, I defined a PetscScalar J[4*n][4*n], and filled in the values for this J matrix by MatSetValues(). What is n? It should not be taking anywhere this much time. How sparse is the matrix? Do you preallocate the nonzero structure? Do you reuse the same matrix for each time step? 245 seconds out of the total 351 seconds in the DAE TS solving part are due to this J matrix computation. For that J matrix, half of them are constants values which doesn't change in each iteration. However, since my J is created inside each IJacobian() call, I couldn't reuse it. If that part of work belongs to redundant computation, I would like to know if there's a way to set up the Jacobian matrix outside of the IJacobian() subroutine, so that I can keep the constant part of values in J for all the iterations but only updates the changing values which depends on X? MatStoreValues() and MatRetrieveValues() but you can only call this after you have assembled the matrix with the correct nonzero structure. So you need to put the constants values in, put zeros in all the locations with non constant values (that are not permeant zeros), call MatAssemblyBegin/End() then call MatStoreValues() then for each computation of the Jacobian you first call MatRetrieveValues() and then put in the non constant values. Then call MatAssemblyBegin/End() Barry Thanks, Shuangshuang -Original Message- From: Jed Brown [mailto:five...@gmail.com] On Behalf Of Jed Brown Sent: Thursday, August 15, 2013 7:27 PM To: Jin, Shuangshuang Cc: petsc-users@mcs.anl.gov Subject: RE: [petsc-users] Performance of PETSc TS solver Jin, Shuangshuang shuangshuang@pnnl.gov writes: Hi, Jed, I followed your suggestion and profiled the IJacobian stage, please see the related profile below: Cool, all of these are pretty inexpensive, so your time is probably in compu
Re: [petsc-users] Performance of PETSc TS solver
90% of the time is spent in your Jacobian evaluation routine which is clearly a bottleneck. On Aug 13, 2013, at 12:34 PM, Jin, Shuangshuang wrote: Hello, Jed and Barry, thanks for your reply. We are solving a power system dynamic simulation problem. We set up the DAE equations and its Jacobian matrix, and would like to use the Trapezoid method to solve it. That's also the reason why we chose TSTHETA. From the PETSc manual, we read that: -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint corresponds to Crank-Nicholson (TSCN). This method can be applied to DAE. For the default Theta=0.5, this is the trapezoid rule (also known as Crank-Nicolson, see TSCN). I haven't heard of ARKIMEX or ROSW before. Are they some external packages or DAE solvers that implement the Trapezoid method? I have also tried the -ksp_type preonly -pc_type lu option you indicated but failed. The PETSC ERROR messages are: No support for this operation for this object type! Matrix format mpiaij does not have a built-in PETSc LU! PETSc does not have a native parallel direct solver. You can use MUMPS (--download-mumps) or superlu_dist (--download_superlu_dist) Attached please see the log_summary for running the TSTHETA with -ts_theta_theta 0.5 and its default ksp solver. Please help me to evaluate the performance and see what's the bottleneck of the slow computation speed. Thanks a lot! Shuangshuang -Original Message- From: Barry Smith [mailto:bsm...@mcs.anl.gov] Sent: Monday, August 12, 2013 6:39 PM To: Jed Brown Cc: Jin, Shuangshuang; petsc-users@mcs.anl.gov Subject: Re: [petsc-users] Performance of PETSc TS solver Also always send the output from running with -log_summary whenever you ask performance questions so we know what kind of performance it is getting. Barry On Aug 12, 2013, at 8:14 PM, Jed Brown jedbr...@mcs.anl.gov wrote: Jin, Shuangshuang shuangshuang@pnnl.gov writes: Hello, PETSc developers, I have a question regarding the performance of PETSc TS solver expecially the TSTHETA. I used it to solve my DAE equations. TSTHETA is not L-stable and not stiffly accurate, so it's not normally something that you'd want to use for a DAE. Make sure you're getting meaningful results and try switching to something like an ARKIMEX or ROSW since those are likely better for your problem. I have recorded the solution times when different numbers of processors are used: 2 processors: 1021 seconds, 4 processors: 587.244 seconds, 8 processors: 421.565 seconds, 16 processors: 355.594 seconds, 32 processors: 322.28 seconds, 64 processors: 382.967 seconds. It seems like with 32 processors, it reaches the best performance. However, 322.28 seconds to solve such DAE equations is too slow than I expected. The number of equations (1152) is quite small, so I'm not surprised there is no further speedup. Can you explain more about your equations? I have the following questions based on the above results: 1. Is this the usual DAE solving time in PETSc to for the problem with this dimension? That depends what your function is. 2. I was told that in TS, by default, ksp uses GMRES, and the preconditioner is ILU(0), is there any other alterative ksp solver or options I should use in the command line to solve the problem much faster? I would use -ksp_type preonly -pc_type lu for such small problems. Is the system dense? 3. Do you have any other suggestion for me to speed up the DAE computation in PETSc? Can you describe what sort of problem you're dealing with, what causes the stiffness in your equations, what accuracy you want, etc. job.out.summary
Re: [petsc-users] how to know the original global index after the partition.
One way to do this is to call ISPartitioningCount() first to get the number of indices on each process and then call ISInvertPermutation() to get the IS that has the original global index. Here's what I've done in one of my code. /* Convert the processor mapping IS to new global numbering */ ierr = ISPartitioningToNumbering(is,is_globalnew);CHKERRQ(ierr); PetscInt *nloc; /* Convert new global numbering to old global numbering */ ierr = PetscMalloc(nps*sizeof(PetscInt),nloc);CHKERRQ(ierr); ierr = ISPartitioningCount(is,nps,nloc);CHKERRQ(ierr); // nps is the number of processors. ierr = ISDestroy(is);CHKERRQ(ierr); ierr = ISInvertPermutation(is_globalnew,nloc[rank],is);CHKERRQ(ierr); The resultant is is the index set in the old global numbering. Shri On Aug 7, 2013, at 5:46 AM, 丁老师 wrote: In the ISPartitioningToNumbering(IS part,IS *is) is define thee index set that defines the global numbers on each part. but how to get the original global index.
[petsc-users] TSStep error
On Mar 27, 2013, at 8:18 PM, Jin, Shuangshuang wrote: Hi, I am using TS to solve a nonlinear DAE problem. I got an error below: [0]PETSC ERROR: - Error Message [0]PETSC ERROR: ! [0]PETSC ERROR: TSStep has failed due to DIVERGED_NONLINEAR_SOLVE, increase -ts_max_snes_failures or make negative to attempt recovery! [0]PETSC ERROR: [0]PETSC ERROR: Petsc Development HG revision: 6e0ddc6e9b6d8a9d8eb4c0ede0105827a6b58dfb HG Date: Mon Mar 11 22:54:30 2013 -0500 [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: ./dynSim on a arch-complex named olympus.local by d3m956 Wed Mar 27 18:12:05 2013 [0]PETSC ERROR: Libraries linked from /pic/projects/ds/petsc-dev/arch-complex/lib [0]PETSC ERROR: Configure run at Tue Mar 12 14:32:37 2013 [0]PETSC ERROR: Configure options --with-scalar-type=complex --with-clanguage=C++ PETSC_ARCH=arch-complex --with-fortran-kernels=generic [0]PETSC ERROR: [0]PETSC ERROR: TSStep() line 2442 in src/ts/interface/ts.c [0]PETSC ERROR: TSSolve() line 2553 in src/ts/interface/ts.c [0]PETSC ERROR: simu() line 420 in unknowndirectory/simulation.C [0]PETSC ERROR: runSimulation() line 83 in unknowndirectory/dynSim.h Run simulation time: 0.0317168 What does this error imply for? The error implies that the nonlinear solve (Newton's method), used at each time step, did not converge. This could be due to various reasons as described here http://www.mcs.anl.gov/petsc/documentation/faq.html#newton My ftime = 0.5, and the timestep is 0.001, and the solution method I?m using is TSBEULER PetscReal ftime=0.5; ierr = TSSetInitialTimeStep(ts, 0.0, 0.001); CHKERRQ(ierr); ierr = TSSetType(ts, TSBEULER); CHKERRQ(ierr); Is there any command option I should use when I run the code? Thanks, Shuangshuang -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130328/81289784/attachment.html
[petsc-users] advice on a linear system solution for a matrix in blocks
PETSc's solver for block matrices PCFieldsplit http://www.mcs.anl.gov/petsc/petsc-3.3/docs/manualpages/PC/PCFIELDSPLIT.html#PCFIELDSPLIT is what you need. More information on fieldsplit is given in section 4.5 of the manual. - Original Message - Dear all, I have a linear system in blocks that I would like to solve. The system is written as [ A C^T ] ? b_1 ? ? ? ? ? ? ? ?= [ C ? ? ? B] ? b_2 where C matrix results from a coupling condition in fluid-structure interaction and rather sparse depending on the size of the coupling surfaces. A and B are also sparse symmetric matrices. I was wondering if I could find a solution to this problem without factorizing the complete operator matrix but only using, somehow, A and B and using some algebraic tricks and the sparsity of C. Any pointers/ideas are appreciated. Best regards, Umut
[petsc-users] PETSc DAE solver
Jin, If you are working on integrating power grid DAE equations then you should use petsc-dev instead as we've added an example for it. petsc-dev/src/ts/examples/tutorials/power_grid/stabiltiy_9bus/ex9bus.c. There are other examples in the power_grid directory. Shri On Feb 20, 2013, at 12:19 PM, Jin, Shuangshuang wrote: Hi, I?m trying to use PETSc to solve a set of first-order differential and algebraic equations such as: x' = f(x,y) 0 = g(x,y) The resources I can find online so far are PETSc Manual: Chapter 6 TS: Scalable ODE and DAE Solvers petsc-3.3-p6/src/ts/examples/tutorials which still seems vague to me at this moment. Can anyone tell me if there?s any more specific references I can refer to as a starting point to tackle this problem? Thanks, Shuangshuang -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130220/f00d35ad/attachment.html
[petsc-users] PETSc DAE solver
http://www.mcs.anl.gov/petsc/developers/index.html On Feb 20, 2013, at 12:39 PM, Jin, Shuangshuang wrote: Shri, thank you! That?s very helpful information. I?m going to obtain the development version of PETSc now. Thanks, Shuangshuang From: petsc-users-bounces at mcs.anl.gov [mailto:petsc-users-bounces at mcs.anl.gov] On Behalf Of Shri Sent: Wednesday, February 20, 2013 10:34 AM To: PETSc users list Subject: Re: [petsc-users] PETSc DAE solver Jin, If you are working on integrating power grid DAE equations then you should use petsc-dev instead as we've added an example for it. petsc-dev/src/ts/examples/tutorials/power_grid/stabiltiy_9bus/ex9bus.c. There are other examples in the power_grid directory. Shri On Feb 20, 2013, at 12:19 PM, Jin, Shuangshuang wrote: Hi, I?m trying to use PETSc to solve a set of first-order differential and algebraic equations such as: x' = f(x,y) 0 = g(x,y) The resources I can find online so far are PETSc Manual: Chapter 6 TS: Scalable ODE and DAE Solvers petsc-3.3-p6/src/ts/examples/tutorials which still seems vague to me at this moment. Can anyone tell me if there?s any more specific references I can refer to as a starting point to tackle this problem? Thanks, Shuangshuang -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130220/a220463e/attachment.html
[petsc-users] PETSc DAE solver
On Feb 20, 2013, at 2:08 PM, Jin, Shuangshuang wrote: Hi, Shri, I was able to run the ex9bus.c code successfully with ?mpirun ?n 1 ex9bus?, but failed with ?mpirun ?n 2 ex9bus? with the following message: [0]PETSC ERROR: - Error Message [0]PETSC ERROR: No support for this operation for this object type! [0]PETSC ERROR: Only for sequential runs! [0]PETSC ERROR: [0]PETSC ERROR: Petsc Development HG revision: fab2442df69c95ef24bf6631a38c7dc028d7ee34 HG Date: Tue Feb 19 16:16:15 2013 -0600 [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: ex9bus on a arch-linux2-c-debug named olympus.local by d3m956 Wed Feb 20 11:36:02 2013 [0]PETSC ERROR: Libraries linked from /pic/projects/ds/petsc-dev/arch-linux2-c-debug/lib [0]PETSC ERROR: Configure run at Wed Feb 20 11:01:15 2013 [0]PETSC ERROR: Configure options [0]PETSC ERROR: [0]PETSC ERROR: main() line 807 in src/ts/examples/tutorials/power_grid/stability_9busex9bus.c Does that mean the parallelized version of DAE solver in PETSc is not available yet? Or I missed some PETSc options to enable the function? This example is written for a single processor only, it cannot be run in parallel. Another question is I also found two other examples in previous installed non-dev PETSc version: petsc-3.3-p6/src/ts/examples/tutorials/ex8.c and ex19.c. These examples are in petsc-dev too. ex8 is actually an ODE but it shows how to use the framework for implicit time-stepping solvers. ex19 is a DAE example. I felt that they are somewhat similar to my problem. Its hard for me to tell what my problem means. What are you trying to solve? Do you think those two are useful examples to me as well? However, it seems like they are uniprocessor examples only too. Yes these are all uniprocessor examples because they are really small. See ex17 for a parallel TS example. Please correct me if I understand it wrong. Thanks, Shuangshuang From: petsc-users-bounces at mcs.anl.gov [mailto:petsc-users-bounces at mcs.anl.gov] On Behalf Of Shri Sent: Wednesday, February 20, 2013 10:49 AM To: PETSc users list Subject: Re: [petsc-users] PETSc DAE solver http://www.mcsanl.gov/petsc/developers/index.html On Feb 20, 2013, at 12:39 PM, Jin, Shuangshuang wrote: Shri, thank you! That?s very helpful information. I?m going to obtain the development version of PETSc now. Thanks, Shuangshuang From: petsc-users-bounces at mcs.anl.gov [mailto:petsc-users-bounces at mcs.anl.gov] On Behalf Of Shri Sent: Wednesday, February 20, 2013 10:34 AM To: PETSc users list Subject: Re: [petsc-users] PETSc DAE solver Jin, If you are working on integrating power grid DAE equations then you should use petsc-dev instead as we've added an example for it. petsc-dev/src/ts/examples/tutorials/power_grid/stabiltiy_9bus/ex9bus.c. There are other examples in the power_grid directory. Shri On Feb 20, 2013, at 12:19 PM, Jin, Shuangshuang wrote: Hi, I?m trying to use PETSc to solve a set of first-order differential and algebraic equations such as: x' = f(x,y) 0 = g(x,y) The resources I can find online so far are PETSc Manual: Chapter 6 TS: Scalable ODE and DAE Solvers petsc-3.3-p6/src/ts/examples/tutorials which still seems vague to me at this moment. Can anyone tell me if there?s any more specific references I can refer to as a starting point to tackle this problem? Thanks, Shuangshuang -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130220/4cda9bc0/attachment-0001.html
[petsc-users] Communication time with threading enabled.
This has something to do with OpenMPI. I cannot reproduce this issue with MPICH. Can you try switching to MPICH (--download-mpich). Shri On Feb 14, 2013, at 6:58 AM, Gon?alo Albuquerque wrote: Dear All, I'm experimenting with PETSc hybrid MPI/OpenMP capabilities and I have run a rather simple test case (2D magnetostatic) using PETSc compiled with both OpenMPI and thread support (both PThreads and OpenMP) on a Ubuntu 12.04 system. I cannot figure out the results obtained when comparing runs made using the same number of MPI processes (2) and specifying either no threads (-threadcomm_type nothread) or 1 OpenMP thread (-threadcomm_type openmp -threadcomm_nthreads 1). I attached the logs of both runs. It seems that the communication time has literally exploded. A grep over the logs gives: No threading: Average time for MPI_Barrier(): 1.38283e-06 Average time for zero size MPI_Send(): 7.03335e-06 1 OpenMP thread: Average time for MPI_Barrier(): 0.00870218 Average time for zero size MPI_Send(): 0.00614798 The same things is occurring when running ksp ex5 (see attached logs). Any ideas as to what I'm missing? Many thanks in advance, Gon?alo nothread.logopenmp_nthreads_1.logex5_nothread.logex5_openmp_nthreads_1.log
[petsc-users] thread model
On Nov 3, 2012, at 8:20 PM, Mohammad Mirzadeh wrote: Hum ... That did not help. I found this -threadcomm_affinities 0,1,2,3,4,5,6,7 but that also did nothing. Is that option any relevant? This option currently only works with pthreads. Another question: the linux box I'm testing on has 2 quad-cores (8 cores/2 sockets). If I run the code with mpirun -np 8, I get about 3X which makes sense to me. However, If I run with either pthread (which seems to use all cores) or openmp (which always defaults to 1 no matter what) As I said before you have to set the core affinity for OpenMP through environment variables. I get the same performance as serial. i) Which example are you running? Please send output of -log_summary. ii) Only the vector and a few matrix operations are threaded currently. There are no threaded preconditioners yet. Does this mean there is something messed up with the hardware and/or how mpi/openmp/pthread is set up? Shri On Fri, Nov 2, 2012 at 8:24 PM, Shri abhyshr at mcs.anl.gov wrote: On Nov 2, 2012, at 9:54 PM, Mohammad Mirzadeh wrote: Hi, To use the new thread model in PETSc, does it suffice to run the code with the following? -threadcomm_type openmp/pthread -threadcomm_nthreads # When I run the code with openmp, only 1 processor/core is active (looking at top). When using pthread, all cores are active. Am I missing something? OpenMP defaults to binding all threads to a single core if the cpu affinity is not specified explicitly. If you are using GNU OpenMP then you can bind the threads to specific cores using the environment variable GOMP_CPU_AFFINITY http://gcc.gnu.org/onlinedocs/libgomp/GOMP_005fCPU_005fAFFINITY.html If you are some other OpenMP implementation then check its manual to see how to set cpu affinity. Shri -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20121103/88ce0677/attachment.html
[petsc-users] thread model
On Nov 2, 2012, at 9:54 PM, Mohammad Mirzadeh wrote: Hi, To use the new thread model in PETSc, does it suffice to run the code with the following? -threadcomm_type openmp/pthread -threadcomm_nthreads # When I run the code with openmp, only 1 processor/core is active (looking at top). When using pthread, all cores are active. Am I missing something? OpenMP defaults to binding all threads to a single core if the cpu affinity is not specified explicitly. If you are using GNU OpenMP then you can bind the threads to specific cores using the environment variable GOMP_CPU_AFFINITY http://gcc.gnu.org/onlinedocs/libgomp/GOMP_005fCPU_005fAFFINITY.html If you are some other OpenMP implementation then check its manual to see how to set cpu affinity. Shri
[petsc-users] Threaded Petsc
Subramanya, Threads can be used in petsc-dev via the thread communicator object threadcomm. PETSc currently has three thread communicators available: i) nothread - This is the default and does not use any threads. ii) pthread - This one uses native lockfree thread pool implementation using POSIX threads iii) OpenMP - OpenMP manages the threadpool. Configuration: Configure with --with-pthreadclasses=1 and/or --with-openmp=1 to use pthreads or OpenMP. **Note: If you are planning to use the OpenMP thread communicator only on MacOS then you'll need to additionally configure with --with-pthreadclasses=1. This is because OpenMP threadprivate variables are not available on MacOS and we use pthread_get/setspecific() calls for thread private variables instead (which are only available when configured with --with-pthreadclasses). Usage: cd src/sys/objects/threadcomm/examples/tutorials make ex1 ./ex3 -threadcomm_type nothread,openmp,pthread -threadcomm_nthreads # of threads Use -h to get the list of available options with threadcomm. Examples: There are several examples in src/sys/threadcomm/examples/tutorials/ .Take a look at these to see how the threadcomm interface is used. We currently have PETSc's vector (seq and mpi) operations threaded but haven't made it open yet. These operations are guarded by the flag PETSC_THREADCOMM_ACTIVE. In order to use the threaded vector operations define PETSC_THREADCOMM_ACTIVE in $PETSC_DIR/$PETSC_ARCH/include/petscconf.h I'll update the docs on the website. Shri On Oct 15, 2012, at 11:59 PM, Barry Smith wrote: Ok, thanks. Shri will update them with the latest info tomorrow. Sorry for the confusion. (they are easier :-)). Barry On Oct 15, 2012, at 11:57 PM, Subramanya G subramanya.g at gmail.com wrote: Hi Barry, These are the instructions that I found. Thanks. http://www.mcs.anl.gov/petsc/petsc-dev/docs/installation.html#threads Subramanya G Sadasiva, Graduate Research Assistant, Hierarchical Design and Characterization Laboratory, School of Mechanical Engineering, Purdue University. The art of structure is where to put the holes Robert Le Ricolais, 1894-1977 On Mon, Oct 15, 2012 at 9:49 PM, Barry Smith bsmith at mcs.anl.gov wrote: On Oct 15, 2012, at 11:35 PM, Subramanya G subramanya.g at gmail.com wrote: Hi, When I try to run the example suggested with -dm_vector_type pthread , I get a unknown vector type error. That documentation is out of date. Please let us know where you saw this suggestion so we can remove it. Barry Unknown vector type: pthread! Here are my configure lines Configure options --download-mpich --download-f-blas-lapack --download-scientificpython --download-fiat --download-generator --download-chaco --download-triangle --with-shared-libraries --download-metis --download-parmetis --download-umfpack --download-hypre --sharedLibraryFlags=-fPIC -shared --with-pthreadclasses I am using the latest Petsc Dev version. Does anybody have any idea what might be wrong? It built without any issues. Do I need to do anything else when I configure. Subramanya G Sadasiva, Graduate Research Assistant, Hierarchical Design and Characterization Laboratory, School of Mechanical Engineering, Purdue University. The art of structure is where to put the holes Robert Le Ricolais, 1894-1977 On Mon, Oct 15, 2012 at 5:44 PM, Barry Smith bsmith at mcs.anl.gov wrote: On Oct 15, 2012, at 7:38 PM, Subramanya G subramanya.g at gmail.com wrote: Hi , I need iterative solvers for some of my code and I was planning to use petsc for it. Owing to the fact that some other parts of my code aren't parallelizable, I can't use petsc in parallel. However, I see that petsc can do threads now.. I have a small question about this . Does this allow me to assemble a matrix in the main thread. and then when solve is called will petsc take care of using multiple threads? Yes Will this give me any speed ups? Yes You need to use petsc-dev for this functionality and let us know if you have difficulties. Barry Thanks/ Subramanya G Sadasiva, Graduate Research Assistant, Hierarchical Design and Characterization Laboratory, School of Mechanical Engineering, Purdue University. The art of structure is where to put the holes Robert Le Ricolais, 1894-1977
[petsc-users] Threaded Petsc
On Oct 16, 2012, at 11:30 AM, Subramanya G wrote: Hi Thanks very much for the info. Does this mean that all the internal vector and matrix operations are multi-threaded? All the vector operations are threaded. All the matrix operations aren't threaded yet (only MatMult and a few others). We are working on threading the matrix operations once we test the vector operations. If I have a sequential code with PETSC_THREADCOMM_ACTIVE set to on will it be faster than a purely sequential code, especially for the KSP routines , which I guess depend on AXPY type operations. Yes, you should expect some speedup for the vector operations and MatMult. Shri Thanks Subramanya G Sadasiva, Graduate Research Assistant, Hierarchical Design and Characterization Laboratory, School of Mechanical Engineering, Purdue University. The art of structure is where to put the holes Robert Le Ricolais, 1894-1977 On Tue, Oct 16, 2012 at 9:10 AM, Shri abhyshr at mcs.anl.gov wrote: Subramanya, Threads can be used in petsc-dev via the thread communicator object threadcomm. PETSc currently has three thread communicators available: i) nothread - This is the default and does not use any threads. ii) pthread - This one uses native lockfree thread pool implementation using POSIX threads iii) OpenMP - OpenMP manages the threadpool. Configuration: Configure with --with-pthreadclasses=1 and/or --with-openmp=1 to use pthreads or OpenMP. **Note: If you are planning to use the OpenMP thread communicator only on MacOS then you'll need to additionally configure with --with-pthreadclasses=1. This is because OpenMP threadprivate variables are not available on MacOS and we use pthread_get/setspecific() calls for thread private variables instead (which are only available when configured with --with-pthreadclasses). Usage: cd src/sys/objects/threadcomm/examples/tutorials make ex1 ./ex3 -threadcomm_type nothread,openmp,pthread -threadcomm_nthreads # of threads Use -h to get the list of available options with threadcomm. Examples: There are several examples in src/sys/threadcomm/examples/tutorials/ .Take a look at these to see how the threadcomm interface is used. We currently have PETSc's vector (seq and mpi) operations threaded but haven't made it open yet. These operations are guarded by the flag PETSC_THREADCOMM_ACTIVE. In order to use the threaded vector operations define PETSC_THREADCOMM_ACTIVE in $PETSC_DIR/$PETSC_ARCH/include/petscconf.h I'll update the docs on the website. Shri On Oct 15, 2012, at 11:59 PM, Barry Smith wrote: Ok, thanks. Shri will update them with the latest info tomorrow. Sorry for the confusion. (they are easier :-)). Barry On Oct 15, 2012, at 11:57 PM, Subramanya G subramanya.g at gmail.com wrote: Hi Barry, These are the instructions that I found. Thanks. http://www.mcs.anl.gov/petsc/petsc-dev/docs/installation.html#threads Subramanya G Sadasiva, Graduate Research Assistant, Hierarchical Design and Characterization Laboratory, School of Mechanical Engineering, Purdue University. The art of structure is where to put the holes Robert Le Ricolais, 1894-1977 On Mon, Oct 15, 2012 at 9:49 PM, Barry Smith bsmith at mcs.anl.gov wrote: On Oct 15, 2012, at 11:35 PM, Subramanya G subramanya.g at gmail.com wrote: Hi, When I try to run the example suggested with -dm_vector_type pthread , I get a unknown vector type error. That documentation is out of date. Please let us know where you saw this suggestion so we can remove it. Barry Unknown vector type: pthread! Here are my configure lines Configure options --download-mpich --download-f-blas-lapack --download-scientificpython --download-fiat --download-generator --download-chaco --download-triangle --with-shared-libraries --download-metis --download-parmetis --download-umfpack --download-hypre --sharedLibraryFlags=-fPIC -shared --with-pthreadclasses I am using the latest Petsc Dev version. Does anybody have any idea what might be wrong? It built without any issues. Do I need to do anything else when I configure. Subramanya G Sadasiva, Graduate Research Assistant, Hierarchical Design and Characterization Laboratory, School of Mechanical Engineering, Purdue University. The art of structure is where to put the holes Robert Le Ricolais, 1894-1977 On Mon, Oct 15, 2012 at 5:44 PM, Barry Smith bsmith at mcs.anl.gov wrote: On Oct 15, 2012, at 7:38 PM, Subramanya G subramanya.g at gmail.com wrote: Hi , I need iterative solvers for some of my code and I was planning to use petsc for it. Owing to the fact that some other parts of my code aren't parallelizable, I can't use petsc in parallel. However, I see that petsc can do threads now.. I have a small question about this . Does this allow me to assemble a matrix in the main thread. and then when solve is called will petsc take care
[petsc-users] Pthread support
Alex, We've taken out the old pthread code since it was experimental and we plan to have a more generic interface that can support different threading models rather than sticking to one. All the thread stuff can now be accessed through the threadcomm interface. If you run your code with --help you'll see the available options available with threadcomm. Shri Shri On Jul 23, 2012, at 7:08 PM, Alexander Goncharov wrote: Hello! I have a question about pthread support in PETSC. I could not find it in the development version. Is it going to be supported in the new release? The reason I tried development version is because in petsc-3.3 compiled with pthreadclasses=1 all jobs would be sent to one core, if I start executables separately. Same for MPI, mpirun -np 4 would show up as 4 processes with 25% CPU usage each, in the top output. Without pthreadclasses all is ok. I would have processes with each utilizing separate core at 100%. Have you come across such a behavior? I ran ex2 from ksp tutorials. Thank you! Alexander.
[petsc-users] Pthread support
On Jul 23, 2012, at 8:00 PM, Alexander Goncharov wrote: Hi Jed, thank you for a quick reply. If I run ex2 with -mat_type seqaijpthread -nthreads 4, everything runs fine. I see 400% utilization. I've taken out all this code from the development version. Are you using petsc-3.3 while running this? The problem persists if I just start two instances of the ex2 from the console without seqaijpthread/seqpthread options. Just two single processor jobs. They show up as two 50% jobs in the top output and run slower. I hope I made my problem more clear. Thank you, Alexander. On Mon, 2012-07-23 at 21:41 -0500, Jed Brown wrote: On Mon, Jul 23, 2012 at 9:08 PM, Alexander Goncharov alexvg77 at gmail.com wrote: Hello! I have a question about pthread support in PETSC. I could not find it in the development version. Is it going to be supported in the new release? The reason I tried development version is because in petsc-3.3 compiled with pthreadclasses=1 all jobs would be sent to one core, if I start executables separately. Same for MPI, mpirun -np 4 would show up as 4 processes with 25% CPU usage each, in the top output. Without pthreadclasses all is ok. I would have processes with each utilizing separate core at 100%. Have you come across such a behavior? Sounds like an affinity problem. When running with multiple MPI processes and threads, you generally have to specify the number of threads manually. The threading code is in the process of being overhauled to a cleaner and more flexible design. It's not ready for use yet, but should be in a couple months.
[petsc-users] Pthread support
The current crude default affinity setting is that for each process the first thread gets assigned to core number 0, second to core number 1, and so on (assuming that the OS numbers the cores sequentially). However, the OS could number the cores differently. For example on Intel Nehalem Architecture, cores on one socket have even numbers (0,2,4,6) and that on the other socket have odd (1,3,5,7). We currently don't have any mechanism by which we can figure out how the OS numbers the cores and set affinity accordingly. Shri On Jul 23, 2012, at 8:33 PM, Barry Smith wrote: On Jul 23, 2012, at 10:22 PM, Jed Brown wrote: On Mon, Jul 23, 2012 at 10:19 PM, Barry Smith bsmith at mcs.anl.gov wrote: Yes, but is PETSc built with MPI? PETSc does not decide on which core to run the program, but MPI might. Well, PETSc is setting affinity if it can when you use threads. Jed, So maybe that is the problem. PETSc is forcing affinity to a hardware number core regardless of whether another PETSc program is running and has also forced that same core? Barry
[petsc-users] Pthread support
I'm at a conference in San Diego so slow in responding to emails. I'll try to get this in tomorrow. Shri On Jul 24, 2012, at 12:16 PM, Jed Brown jedbrown at mcs.anl.gov wrote: On Tue, Jul 24, 2012 at 2:13 PM, Alexander Goncharov alexvg77 at gmail.com wrote: I actually tried taskset utility, and it did not help. As explained in the previous email, affinity is naively done by calling CPU_ZERO and then adding in one core per thread. Whenever my suggestion gets implemented, threads will inherit from taskset. For now, you can use -threadcomm_affinities ..., but that's low-level enough to be painful. -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120724/b0923b92/attachment-0001.html
[petsc-users] How do you write/read to the .info file
You can also use the option -viewer_binary_skip_info to have PETSc not create the .info file. - Original Message - On Tue, May 22, 2012 at 8:45 PM, Andrew Spott andrew.spott at gmail.com wrote: When writing out a Mat in binary form, the .info file is created, but it usually empty, how do you read/write to it? They are written and read automatically by the PETSc binary viewers. You shouldn't need to look at them. -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120522/117f0c97/attachment.html
[petsc-users] Modifying the structure of a matrix.
On May 14, 2012, at 11:07 PM, Andrew Spott wrote: I'm attempting to do a convergence study on my basis, see how large my basis need to be to make my results accurate, does that help? You could use create a matrix with the set of all basis vectors and then zero out the rows columns of the matrix using MatZeroRowsColumns. http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatZeroRowsColumns.html I think this should be better than each time creating a new matrix when the basis changes. -Andrew On May 14, 2012, at 9:55 PM, Shri wrote: What are you trying to do by removing the rows/columns of the matrix? Are there any variables that you want to treat as constants and hence eliminating the rows and columns? Shri On May 14, 2012, at 10:51 PM, Andrew Spott wrote: Is there a computational cost to using a sparse matrix that is much bigger than it should be? for example, say you have a matrix that is n x n, but you only store anything in the first m rows and columns (m n). Also, would you then have to redo the matrix distribution among nodes by hand, or will PETSC do it for you? Thanks for the quick reply -Andrew On May 14, 2012, at 9:38 PM, Jed Brown wrote: On Mon, May 14, 2012 at 9:33 PM, Andrew Spott andrew.spott at gmail.com wrote: What is the best way to selectively add/remove columns/rows of a Mat? Is there a way to do it in place, or does it require essentially moving the values to a new matrix? If you only need to add/remove a small number of values, it's best to allocate for all that might be present and just store explicit zeros. To add or remove rows, you are also changing the vector size that can be multiplied against, and the data structure must be rebuilt. This is why it's frequently preferable to use MatZeroRows (or a variant) instead of removing them from the data structure. -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120514/b3548b79/attachment-0001.htm
[petsc-users] Modifying the structure of a matrix.
By default, PETSc uses a row distribution with the first m rows on processor 0, next m rows on processor 1, and so on. If the order is not important, then you could remove rows/columns from each processor instead of removing the last n-m rows on the last processor. Thus, the load balance would still be even. Shri - Original Message - I'm solving the schrodinger equation in a basis state method and looking at adding the azimuthal quantum numbers (the m's) to my basis so I can look at circularly polarized light. However, when I do this, I'll need to do some sort of convergence study to make sure I add enough of them. The way my code is structured, it will probably be easier to just remove rows and columns from a bigger matrix, instead of adding them to a smaller matrix. However, depending on the way I structure the matrix, I could end up removing all the values (or a significant portion of them) from a processor when I do that. Speaking more on that, is the PETSC_DECIDE way of finding the local distribution smart in any way? or does it just assume an even distribution of values? (I assume that it assumes a even distribution of values before the assembly, but does it redistribute during assembly?) Thanks, -Andrew On May 15, 2012, at 6:40 AM, Jed Brown wrote: On Mon, May 14, 2012 at 10:39 PM, Andrew Spott andrew.spott at gmail.com wrote: That is what I figure. I'm curious though if you need to manually determine the local row distribution after you do that. (for example, say you completely remove all the values from the local range of one processor? that processor wouldn't be utilized unless you redistribute the matrix) What sizes and method are we talking about? Usually additional (compact) basis functions only make sense to add to one of a small number of processes. -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120515/1d82f08e/attachment.htm
[petsc-users] Modifying the structure of a matrix.
What are you trying to do by removing the rows/columns of the matrix? Are there any variables that you want to treat as constants and hence eliminating the rows and columns? Shri On May 14, 2012, at 10:51 PM, Andrew Spott wrote: Is there a computational cost to using a sparse matrix that is much bigger than it should be? for example, say you have a matrix that is n x n, but you only store anything in the first m rows and columns (m n). Also, would you then have to redo the matrix distribution among nodes by hand, or will PETSC do it for you? Thanks for the quick reply -Andrew On May 14, 2012, at 9:38 PM, Jed Brown wrote: On Mon, May 14, 2012 at 9:33 PM, Andrew Spott andrew.spott at gmail.com wrote: What is the best way to selectively add/remove columns/rows of a Mat? Is there a way to do it in place, or does it require essentially moving the values to a new matrix? If you only need to add/remove a small number of values, it's best to allocate for all that might be present and just store explicit zeros. To add or remove rows, you are also changing the vector size that can be multiplied against, and the data structure must be rebuilt. This is why it's frequently preferable to use MatZeroRows (or a variant) instead of removing them from the data structure. -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120514/ab63e191/attachment.htm
[petsc-users] Complex Number
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscOptionsGetScalar.html - Original Message - Dear users n devlopers, My Petsc routine contans a complex number, as a parameter of my problem. I need to provide it at the run time options, since I intend to test a variety of parameter. In this regard, I read that PetscScalar is the way to declare a complex number. If I need merely a complex number at run time option to compile a compact routine, any details from any examples of libraries ? Thanks in anticipation. Abdul, -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120416/51a5e59b/attachment.htm
[petsc-users] Solving a Linear System with two different Solvers (one each thread)..?
Pedro, Taking in to account that pthread was recenltly added to petsc, is it possible to code a program that solve the same linear system with two different solvers, one on each thread?. I guess that can be done calling different ksp object, one on each thread, but I don't have clear, first if this is possible and second how to do that. No, this is not possible currently and may not be available in the future too. I don't think there may be any performance benefit by solving the same linear system using different KSPs on different threads. What are you trying to do using pthreads in PETSc? Shri - Original Message - Hello Petsc Team, I has some experience with PETSc + MPI, but neither with PETSc+Pthreads nor Pthreads alone. Taking in to account that pthread was recenltly added to petsc, is it possible to code a program that solve the same linear system with two different solvers, one on each thread?. I guess that can be done calling different ksp object, one on each thread, but I don't have clear, first if this is possible and second how to do that. Sorry if my question is meaningless, I just started to study pthreads. I will really appreciated any advice. Thanks in advance. Best Regards -- Pedro Torres -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120319/35dbe609/attachment.htm
[petsc-users] Petsc Performance Help
As Jed points out, most of the time (~98%) is spent in your code rather than in the PETSc solver. Clearly, you need to profile your own source code and figure out where the most time is spent. You can do this using the tools for user code profiling (PetscLogEvent and PetscLogStage) in PETSc. Read Sections 11.2 and 11.3 of the User's Manual for more information on PetscLogEvent and PetscLogStage and also look at examples such as http://www.mcs.anl.gov/petsc/petsc-current/src/vec/vec/examples/tutorials/ex5.c.html http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex2.c.html http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex46.c.html Shri - Original Message - On Fri, Mar 16, 2012 at 14:52, Nan.Jia at dartmouth.edu wrote: Dear PETSc Group, I am tuning the efficiency of my PETSc code for a while, but get very little progress. So can anyone help me to analysis the log? Any suggestions will be appreciated. My problem is time dependent. At every time step, two about 6000 by 6000 sparse matrices need to be solved, which come from a Poisson equation. I use both sequential and parallel AIJ format to store matrices, but the performances are both not very good. 1. You need to heed the huge warning ## # # # WARNING!!! # # # # This code was compiled with a debugging option, # # To get timing results run config/configure.py # # using --with-debugging=no, the performance will # # be generally two or three times faster. # # # ## 2. You only spend 18 of 880 seconds in the solver. What do you want? 3. The problem is too small to get significant parallel speedup. http://www.mcs.anl.gov/petsc/documentation/faq.html#slowerparallel -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120316/6fc6fe8d/attachment.htm
[petsc-users] Sparse matrix partitioning in PETSc
Does your matrix come from an unstructured grid? If so, then you'll need to use a partitioning package such as permits. Read section 3.5 of the manual http://www.mcs.anl.gov/petsc/petsc-current/docs/manual.pdf and also see the example http://www.mcs.anl.gov/petsc/petsc-current/src/mat/examples/tutorials/ex15.c.html - Original Message - I want to load a sparse matrix from a file (in PETSc binary format). The matrix will be in MPIAIJ format. I want it is a balanced load: each processor will have nearly equal number of nonzeros. I don't come up with a sequence of PETSc calls to do that. On Sat, Jan 14, 2012 at 11:41 AM, Jed Brown jedbrown at mcs.anl.gov wrote: On Sat, Jan 14, 2012 at 12:38, Junchao Zhang junchao.zhang at gmail.com wrote: Does PETSc provide convenient functions to compute this layout( i.e.,# rows on each processor), or I have to do it myself? I browsed PETSc document and did not find them. MatSetSizes() lets you specify local and/or global sizes. If you specify both, it checks that they are compatible. The implementation uses http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscSplitOwnership.html If these are not suitable for your purposes, can you be more specific about what you would like to do? -- Junchao Zhang -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120114/c7212e84/attachment.htm
[petsc-users] Sparse matrix partitioning in PETSc
MatLoad() does not distribute the rows based on the number of non zeros. You'll need to first compute the number of rows on each process that gives you equal/nearly equal number of non zeros and then call MatSetSizes() followed by MatLoad(). - Original Message - Hello, When I use MatLoad() to load a sparse matrix, I want each process to have equal number of nonzeros, instead of equal number of rows. How could I achieve that in PETSc? Thanks! -- Junchao Zhang -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120113/49aa5819/attachment.htm
[petsc-users] Return value of PCFactorGetMatrix
MatGetFactor() only sets up the matrix for factorization but does not do any factorization. You need to call MatCholeskyFactorSymbolic() http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCholeskyFactorSymbolic.html and MatCholeskyFactorNumeric() http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCholeskyFactorNumeric.html#MatCholeskyFactorNumeric after MatGetFactor(). Shri - Original Message - Dear list, following advice from this list and petsc documentation (especially http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex52.c.html ), i have the following code: //- ierr = KSPCreate( coml,kspBA);CHKERRQ(ierr); ierr = KSPSetOperators( kspBA, Kt, Kt, DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); ierr = KSPSetType(kspBA,KSPPREONLY);CHKERRQ(ierr); ierr = KSPGetPC(kspBA, pcba); ierr = PCSetType(pcba, PCCHOLESKY);CHKERRQ(ierr); ierr = PCFactorSetMatSolverPackage(pcba, MATSOLVERPETSC);CHKERRQ(ierr); ierr = PCFactorSetUpMatSolverPackage(pcba);CHKERRQ(ierr); ierr = MatGetFactor(Kt,petsc,MAT_FACTOR_CHOLESKY,DiagM);CHKERRQ(ierr); ierr = PCFactorGetMatrix(pcba, DiagM);CHKERRQ(ierr); ierr = KSPSetUp(kspBA);CHKERRQ(ierr); //ierr = MatView(DiagM,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); //ierr = MatView(DiagM,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); ierr = VecDuplicate(fin,diag);CHKERRQ(ierr); ierr = MatGetDiagonal(DiagM,diag);CHKERRQ(ierr); //--- where Kt is a stiffness matrix. I was hoping that PCFactorGetMatrix() makes DiagM become D from a LDL^T factorization, but to me it seems as if this is not true. Can anyone clarify? Regards, Uwe
[petsc-users] Conditional Constraints
Jon, How could I implement this constraint inside the PETSc solution loop? You can set bounds on variables directly through TSVISetVariableBounds() http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-dev/docs/manualpages/TS/TSVISetVariableBounds.html Is there any advantage to using SNESVI when my update step is explicit? Or is there a feature I could use to add a Vmag_new computation after each SNES step, and have this taken into account for the SNES convergence criteria? Unfortunately, we do not have some routine like SNESPostIteration(). However, you do this in your nonlinear function evaluation routine by having two variables for SNES iteration numbers i (current iteration) and i_old (previous iteration) i) Get the current SNES iteration number i ii) If i i_old then update Vmag_new and i_old. - Original Message - Hi Shri, I have continued to have little success with my SNESVI formulation of the current constraint. However, working outside of PETSc, I believe I have found an explicit formulation that works well: 1) Set Vmag = Vmag_max and use PETSc to solve for the electric potential distribution (V_ijk is the only degree of freedom). 2) For each electrode, let Vmag_new = Imag_max * Vmag_old / Imag_old. 3) Repeat until Imag is sufficiently close to Imag_max, or until Vmag = Vmag_max with Imag Imag_max. This process converges in about 5 iterations, which is great, but I would like to incorporate the thermal problem again (including TS) now that this electrical problem is solved. How could I implement this constraint inside the PETSc solution loop? Is there any advantage to using SNESVI when my update step is explicit? Or is there a feature I could use to add a Vmag_new computation after each SNES step, and have this taken into account for the SNES convergence criteria? Thanks again for your continued patience and interest in helping me with this problem. Jon On 2011-09-15, at 12:07 PM, Shri wrote: Jon, Additionally, the potential distribution does not look correct as it did before (even without SNES reporting convergence)---instead, the magnitude of the distribution looks like the real part of the correct distribution, if that makes sense. Hmm, in the SNESVI formulation using complex numbers, we have the check real(xlower) = real(x) = real(xupper). i am not sure whether this is the cause for such a behavior. I've had conversation with the other PETSc developers to use absolute values for comparison instead of real values, i.e. abs(xl) = abs(x) = abs(xupper). By doing this, you can essentially have a function for Ve and its bounds being 0 and Vmag_max. However, the other developers think that absolute values being nonlinear in nature might produce some odd behavior in the solution process, so i've havent' support for it yet. This asymmetry may be due to the fact that I now only record the Vmag and Imag variables for the top block of each electrode, whereas before I recorded Ve for all three blocks. Yes, do you want the magnitude and angle for all the three blocks to be same? If so, then probably you'll have to add another constraint. If a block (with coordinates i, j, k) is not occupied by an electrode , C_e is set to zero in the expression for Vb_ijk and all function and Jacobian entries for Vmag and Imag are either not set or set to zero. Do you already know the Vmag and Imag values are these blocks? Are they zero/non-zero? I suggest the following for these blocks xl_ijk_DOF_VMAG = Vmag0_ijk xu_ijk_DOF_VMAG = Vmag0_ijk xl_ijk_DOF_IMAG = Imag0_ijk xu_ijk_DOF_IMAG = Imag0_ijk where Vmag0_ijk and Imag0_ijk is the initial voltage and current magnitudes at these blocks. Basically, the idea here is to keep these magnitudes constant and not have the functions for these variables be a part of the solution process. Shri - Original Message - Hi Shri, Thanks for following up with me. I have made some progress, although my problem is not yet solved. Inspired by your recent suggestions, I realized Ve was redundant as a degree of freedom and did not need to be computed separately from Vmag. I do not want the phase of Ve to be modified by the solution process anyways, so whenever I need the complex value of Ve, I can simply compose it from the Vmag and the constant phase I have recorded for a given electrode. My new function and Jacobian formulations are below if you are interested in seeing them. After I made this change, I began to see the desired behaviour where Vmag decreases when Imag is too high. However, it seems to converge on some magnitude that is maybe 200% too high, and the SNES reports DIVERGED_LINE_SEARCH. Additionally, the potential distribution does not look correct as it did before (even without SNES reporting convergence)---instead, the magnitude of the distribution looks like the real part of the correct distribution, if that makes
[petsc-users] Conditional Constraints
Would f(Ve_ijk) = Imag - Imag_max = 0 simply be equivalent to saying, solve the real system of equations Imag - Imag_max = 0 and 0 = 0 (the imaginary part)? Yes. Regardless, I am not sure what a correct choice of f(Ve_ijk) would be Maybe, having the equation for current block iE_block_ijk = aE_ijk * (Ve - Vb_ijk) for f(Ve_ijk)? where iE_block_ijk = Imag_max + j*0 ? Perhaps, you can have a different and better expression for f(Ve_ijk) based on some related physics. - Original Message - Shri, Note here that (a) f(Ve_ijk) is a complex function while Imag and Imag_max are real. So the above equation does not hold. Would f(Ve_ijk) = Imag - Imag_max = 0 simply be equivalent to saying, solve the real system of equations Imag - Imag_max = 0 and 0 = 0 (the imaginary part)? Regardless, I am not sure what a correct choice of f(Ve_ijk) would be. (b) In your previous email, you had bounds on Imag, 0 Imag Imag_max. Now you want Ve_ijk to increase if Imag Imag_max and decrease if Imag Imag_max. Are you trying to find the value of Ve at which the current is maximum, i.e. Imag = Imag_max, by doing this? Yes, that is essentially what I am trying to find. Ve determines Vb, and the difference between the two determines Imag. Vmag (magnitude of Ve) and Imag have the restrictions 0 Vmag Vmag_max and 0 Imag Imag_max. I want to find either Ve for Imag = Imag_max or Imag for Vmag = Vmag_max. In my test scenario, Vmag_max happens to be quite high so Imag_max will always be reached first, but in general either situation is possible. Thank you again, Jon On 2011-09-08, at 8:34 PM, Shri wrote: Jon, I do not want Ve to be constant, however; I would like it to either decrease in response to Imag Imag_max or increase in response to Imag Imag_max, providing the condition Vmag = Vmag_max remains true. Should I then have something like f(Ve_ijk) = Imag - Imag_max? Note here that (a) f(Ve_ijk) is a complex function while Imag and Imag_max are real. So the above equation does not hold. (b) In your previous email, you had bounds on Imag, 0 Imag Imag_max. Now you want Ve_ijk to increase if Imag Imag_max and decrease if Imag Imag_max. Are you trying to find the value of Ve at which the current is maximum, i.e. Imag = Imag_max, by doing this? - Original Message - Shri, Thanks very much for your response. In fact, -snes_converged_reason revealed that the solution was not converging, as you suggested: Nonlinear solve did not converge due to DIVERGED_LINE_SEARCH. I had not realized this earlier since the voltage distribution in the solution appeared to be correct, even though the SNES_VI constraints were not satisfied. Given that and your explanation, I am sure something is wrong with my choices of f. I do not want Ve to be constant, however; I would like it to either decrease in response to Imag Imag_max or increase in response to Imag Imag_max, providing the condition Vmag = Vmag_max remains true. Should I then have something like f(Ve_ijk) = Imag - Imag_max? Note that you have One more thing from your response: If a block (with coordinates i, j, k) is not occupied by an electrode , C_e is set to zero in the expression for Vb_ijk and all function and Jacobian entries for Ve, Vmag, and Imag are either not set or set to zero. Does this mean f(Ve_ijk) = f(Vmag_ijk) = f(Imag_ijk) = 0? If so, then Ve_ijk, Vmag_ijk, and Imag_ijk would retain their initial values for all time steps. In my present scheme, f(Ve_ijk) = f(Vmag_ijk) = 0 and f(Imag_ijk) is non-zero for electrode blocks; all three are zero for non-electrode blocks. Depending on the choice of initial conditions, I can make f(Vmag_ijk) non-zero as well (until Ve changes to match Ve_max). Thanks again, Jon On 2011-09-08, at 3:42 PM, Shri wrote: Jon, for Ve_ijk, f = 0.0 + 0.0*i, and all corresponding Jacobian entries are zero as well. Ve_ijk is set to the maximum voltage for each electrode when the initial solution is - set. Since the magnitude of Ve is constrained rather than Ve itself, I believe this is correct, but could this be why Ve does not change as it should when the problem is - being solved? This f(Ve_ijk) seems very odd to me since normally you would always have f(some/all variables) = 0 as the function. If f=0.0+0.0*i for all steps, then you don't need a seperate variable Ve_ijk, you can treat it as a constant since Ve_ijk will not change its value during the solutions. Moroever, the Jacobian rows for df(Ve_ijk)/d(all variables) would be all zeros which implies a singular matrix. I am not sure how the solution converged for this case. Were you using -pc_type lu with some shift or some external direct solver such as superlu? What do you get for -snes_monitor -snes_converged_reason? Ve_ijk is set to the maximum voltage for each electrode when
[petsc-users] Iterative solver in Petsc on a multiprocessor computer
Amrit, These plots are just reiterating the problem you've stated previously but provide no clue as to what's the issue here. I suggest the following : Try running your code with a direct solver, mpiexec -n 6 test3 -ksp_type preonly -pc_factor_mat_solver_package mumps, and see if you are getting the expected results on 6 processors. If you do get the expected results using a direct solver then, as Jed points out, the iterative solver options need to be tuned as Jed points. Note that while a direct solver may be able to solve your problem it is not scalable hence iterative solvers are prefered. If the results are not correct then perhaps something is amiss with your code. The default 'parallel' preconditioner in PETSc is the block-jacobi preconditioner with an ILU(0) factorization which may not be a suitable preconditioner on more number of processors. Try mpiexec -n 6 test3 -ksp_type gmres -pc_type bjacobi -sub_pc_type lu This does an LU factorization on each block instead of ILU(0) Perhaps, you should also attach the code, and the command line options. Shri - Original Message - Please find the attached plots where I've compared results obtained from 2 processors (2.pdf) numerical simulation in PETSc with analytical result and they agree fine, but for 6 processors (6.pdf) simulation, the numerical result is a complete garbage. My major confusion is this : why does the default GMRES *works* for 2 processors and not for 6 processors? GMRES algorithm is definitely not the problem here because otherwise it would not give me correct result when I run my simulation on 2 processors. Please note that I *do not* change anything in my code other than the number of processors with this command at run time : mpiexec -n 2 test3 Also, I make sure that I write the output to a separate file. It is very hard for me to believe when you say : (a) the problem you are building is not actually the same or the results are being misinterpreted or It is the same. Like I said above, I am only changing the number of processors at run time. (b) the default methods are not working well and needs some algorithmic adjustments for your problem. if this were an issue, then why does it work for two processors? The problem I am simulating is a vector wave equation in two dimensions (or in other words, Maxwell's equations). The coefficient matrix is assembled using finite difference frequency domain method on a Yee grid. I will appreciate any help in the regard. Thanks. -Amrit Date: Thu, 8 Sep 2011 20:07:04 +0200 From: jedbr...@mcs.anl.gov To: petsc-users at mcs.anl.gov Subject: Re: [petsc-users] Fwd: nonzero prescribed boundary condition On Thu, Sep 8, 2011 at 19:57, amrit poudel amrit_pou at hotmail.com wrote: After running my simulation multiple times on a multiprocessor computer I've just verified that using iterative solver (default gmres) in PETSc to solve a linear system of equations ( Cx=b) with more than 2 processors setting ALWAYS lead to erroneous result. Running identical code with identical setting except for the number of processors ( set this to 2) ALWAYS gives me correct result . You have to explain what erroneous result means here. I am really not sure what is the point behind including iterative solvers if they result into erroneous result on a multiprocessor computer. The result I get from multiprocessor computer is a complete garbage, so I am really not talking about small percentage of error here. Also, if somebody could enlighten why the iterative solvers are error prone on multiprocessors that will be highly appreciated. Well let's not jump to conclusions. Iterative solvers can fail, as can direct solvers, but it's more common that (a) the problem you are building is not actually the same or the results are being misinterpreted or (b) the default methods are not working well and needs some algorithmic adjustments for your problem. Please explain what kind of problem you are solving, how you are going about it, and what symptoms you have observed. -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110908/8eb2ff5c/attachment.htm
[petsc-users] Conditional Constraints
- Original Message - On Sun, Jul 31, 2011 at 22:11, Vijay S. Mahadevan vijay.m at gmail.com wrote: I saw from the documentation that the GL methods were capable of handling DAE systems. But the intention of my question was about the usage of the newton solver with functional constraints (with VI) for a DAE since in a sense, the DAE systems are a composite of an ODE and a constraint. So I was curious whether we could use temporal methods that were not intended for DAEs but for stiff ODEs only and use the VI solver to augment/restrict the solution with constraints. The reason for this is due to the fact that in the creation of a stiff solver, there are additional coefficient constraints that need to be satisfied to be used for a DAE but this is less restrictive for solving ODEs. Do you want equality or inequality constraints? Equality constraints are naturally handled by writing a DAE. Inequality constraints would be handled by VI. VI is not currently available through TS, but my understanding is that Shri is looking into this and it will eventually be available there too. I have already added TSVISetVariableBounds() by which lower and upper bounds can be set for variables directly through the TS interface. I am currently working on getting an example ready for TS with VI solver. Hopefully, i will add it by this week. Thanks, Shri -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110731/c23d6201/attachment.htm
[petsc-users] Conditional Constraints
Jon,
I am trying to compile petsc with complex data type right now and will try
to test the VI solver with complex variables within the next few hours.
My understanding is that comparison operators (such as and ) are
not defined for complex numbers, so would it be correct to say that
anyone using complex scalars within PETSc will encounter this problem
with SNESVI? In other words, SNESVI inequality constraints can only be
applied to real variables or real-valued functions of complex
variables?
For complex scalars, we only compare the real part of the complex variable and
bound,i.e.,
real(xl) = real(x) = real(x_u). For your application, this real part
comparison is sufficient
since the magnitude is saved in the real part.
The residuals are set as follows, and the Jacobian entries are set
correspondingly:
f[V_block] = {function of adjacent V_blocks, V_electrode}
f[V_electrode] = 0.0 + 0.0i
f[magnitude(V_electrode)] = 0.0 + 0.0i
f[magnitude(I_electrode)] = 0.0 + 0.0i
Can you give an example of f for all dofs maybe just for one grid point with
the following dof names
Vb - V_block
Ve - V_electrode
Vmag - magnitude(V_electrode)
Imag - magnitude(I_electrode)
Since I set all the residuals for both magnitude variables to zero,
and each entry in the Jacobian matrix is a partial derivative of a
residual, I believe the Jacobian entries for both magnitudes would all
be zero as well. However, your point made me realize I have no
Jacobian entry for the partial derivative of f[V_block] w.r.t.
V_electrode. Is that needed?
Since f[V_block] is a function of V_electrode, you would need to set the
partial derivatives of f[V_block] w.r.t V_electrode
for correct jacobian evaluation and for having a good newton convergence.
How would one enter it with
MatSetValuesStencil?
ex28.c,
http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-dev/src/ksp/ksp/examples/tutorials/ex28.c.html,
is a petsc example which shows how to use MatSetValuesStencil with multiple
degrees of freedom.
See the routine ComputeMatrix in the example.
My working understanding was that each degree of
freedom conceptually has its own residual and Jacobian, even though
all residuals are contained in one Vec and all Jacobians contained in
one Mat.
Please bare with us as we are trying to figure out whether the issue is with
our VI implementation or something in
your code is amiss. I appreciate your patience on this matter.
Btw, one quick question, are all the problem formulations in your field done
using complex variables or by decomposing the complex term into real/imaginary
or magnitude/angle? My background is in electrical power grid and we mostly use
either real/imaginary part or magnitude/angle as the variables.
Shri
- Original Message -
Thanks for your reply, Shri. Responses below:
On 2011-07-28, at 1:03 PM, Shri wrote:
- Original Message -
Hi Barry, Shri,
I reworked my implementation to include only the time-independent
electrical problem, and to make use of SNESVISetVariableBounds().
However, it does not work quite as I expected: Constraints on
variables that are functions of other variables have no effect on
the
variables by which the functions are defined.
My constraints are set on the four degrees of freedom as follows:
-SNES_VI_INF = V_block = SNES_VI_INF
-SNES_VI_INF = V_electrode = SNES_VI_INF
0.0 + 0.0i = magnitude(V_electrode) = magnitude(V_electrode_max)
0.0 + 0.0i = magnitude(I_electrode) = magnitude(I_electrode_max)
Note that while V_block and V_electrode and complex quantities,the
magnitude terms are actually real numbers. What function did you use
to compute the magnitude? abs()? sqrt of product of complex number
with its conjugate?
I used the 0.0 + 0.0i there to emphasize that I built PETSc with
complex scalars, so all the scalars are stored as complex even if the
imaginary components are always zero. I have been using std::abs()
from the C++ standard library complex template class
(std::complexdouble). I believe the values of type 'double' returned
by std::abs() are implicitly cast to the real part of
std::complexdouble or PetscScalar, though I usually try to cast
explicitly.
This seems like an application where the variables are of type
complex but the constraints need to be set for functions of the
variables which are actually real numbers. I am not sure whether
this can be supported. Maybe Barry or others can shed more
light on this.
My understanding is that comparison operators (such as and ) are
not defined for complex numbers, so would it be correct to say that
anyone using complex scalars within PETSc will encounter this problem
with SNESVI? In other words, SNESVI inequality constraints can only be
applied to real variables or real-valued functions of complex
variables?
The residuals are set as follows, and the Jacobian entries are set
correspondingly:
f[V_block
[petsc-users] Conditional Constraints
- Original Message -
Hi Barry, Shri,
I reworked my implementation to include only the time-independent
electrical problem, and to make use of SNESVISetVariableBounds().
However, it does not work quite as I expected: Constraints on
variables that are functions of other variables have no effect on the
variables by which the functions are defined.
My constraints are set on the four degrees of freedom as follows:
-SNES_VI_INF = V_block = SNES_VI_INF
-SNES_VI_INF = V_electrode = SNES_VI_INF
0.0 + 0.0i = magnitude(V_electrode) = magnitude(V_electrode_max)
0.0 + 0.0i = magnitude(I_electrode) = magnitude(I_electrode_max)
Note that while V_block and V_electrode and complex quantities,the magnitude
terms are actually real numbers. What function did you use to compute the
magnitude? abs()? sqrt of product of complex number with its conjugate?
This seems like an application where the variables are of type complex but the
constraints need to be set for functions of the variables which are actually
real numbers. I am not sure whether this can be supported. Maybe Barry or
others can shed more
light on this.
The residuals are set as follows, and the Jacobian entries are set
correspondingly:
f[V_block] = {function of adjacent V_blocks, V_electrode}
f[V_electrode] = 0.0 + 0.0i
f[magnitude(V_electrode)] = 0.0 + 0.0i
f[magnitude(I_electrode)] = 0.0 + 0.0i
What equations do you use for the partial derivatives of the magnitude
functions w.r.t. the actual complex variables (V_block and V_electrode) in the
jacobian evaluation?
The first two potentials, V_block and V_electrode, are independent,
while magnitude(V_electrode) is a function of V_electrode only and the
electrode current magnitude(I_electrode) is a function of V_block and
V_electrode. Note that V_electrode acts as a source term in the finite
difference formulation for V_block. While all of these variables are
complex, the latter two always have zero imaginary parts. I suppose I
am assuming that PETSc knows to compare only the real parts if the
imaginary parts are zero. (Is this a bad assumption?)
When solving with PETSc (using -snes_vi_type rs)
Did you also use -snes_type vi?
V_electrode stays at
its maximum, as set in the initial guess, and the V_block distribution
falls properly from that. However, measurements of
magnitude(I_electrode) are way above maximum and do not move downwards
with successive SNES iterations. I can also set the constraint on
magnitude(V_electrode) to less than maximum and it does not affect the
value of V_electrode. How can I tell PETSc to change V_electrode when
the magnitude(V_electrode) or magnitude(I_electrode) constraints are
not met?
Send the code for your application with instructions on how to run it and
we'll try to figure out what needs to be
done.
Shri
Thanks again for your help.
Jon
On 2011-07-25, at 8:58 PM, Barry Smith wrote:
On Jul 25, 2011, at 5:50 PM, Jonathan Backs wrote:
Hi Shri,
Thanks for your message and all the helpful tips. If the
TSVISetVariableBounds() functions are available now, I would like
to try them as well. Is the interface essentially the same as with
SNESVISetVariableBounds()? I will get back to you and Barry when I
have had a chance to modify my application.
For my problem, I believe it makes sense to have bounds on one of
the variables as well as one function of the variables. The two
relevant degrees of freedom are the block voltage (one for each
finite difference block) and the electrode voltage (one for each
electrode, which may be present in multiple blocks). The electrode
voltage should keep a constant phase while the magnitude is
constrained between zero and some maximum. The block voltages near
the electrodes depend on the electrode voltages as well as the
neighbouring block voltages. The electrode current for a given
electrode is a function of its electrode voltage and several nearby
block voltages, and should be constrained in magnitude between zero
and some maximum (and the phase unconstrained). Would I need to use
SNESVISetComputeVariableBounds since the electrode current is a
function of the other variables? Would any other provisions need to
be made for the block voltages, since they depend on the electrode
voltages?
You cannot directly make a bound on some function of other
variables. Instead you need to introduce another variable that is
defined to be equal to that function and put the bound on that new
variable.
Barry
Thank you again,
Jon
On 2011-07-25, at 11:58 AM, Shri wrote:
- Original Message -
On Jul 22, 2011, at 4:16 PM, Jonathan Backs wrote:
Barry,
Thank you so much for your response. Lucky, indeed! I look
forward
to trying out these new features.
I neglected to mention in my original post that my electrical
problem is part of a DAE, which includes a time-dependent
[petsc-users] Conditional Constraints
- Original Message - Barry/Shri, Is the new VI SNES type in general terms a constrained nonlinear solver ? Yes, where the constraints are on the variables. Given a function f(x) with upper and lower bounds on x, the VI solver tries to compute an x* such that for each x_i one of the following conditions is satisfied f(x_i) = 0 with xl_i x_i xu_i f(x_i) 0 with x_i = xl_i f(x_i) 0 with x_i = xu_i In other words, given a certain constraint as a function of the variables of interest, does it find the optimal solution that both minimizes the residual and satisfy the constraint ? I would say that the VI solver is not an optimizer but rather a constrained nonlinear solver. The only constrained optimization it can do is if you provide the gradient as f(x) with lower/upper bounds on x. I have been looking at L-BFGS constrained optimization methods and from what I understand, the variational inequality solver seems to fall under a similar but much broader category. I suggest you take a look at the TAO package if you are interested in optimization related problems. http://www.mcs.anl.gov/research/projects/tao/ I would much appreciate it if you can point me to some papers related to the method. http://www.mcs.anl.gov/~tmunson/papers/semismooth.pdf And if I misunderstood the domain of applicability, please do feel free to correct me. Thanks, Vijay On Mon, Jul 25, 2011 at 9:58 PM, Barry Smith bsmith at mcs.anl.gov wrote: On Jul 25, 2011, at 5:50 PM, Jonathan Backs wrote: Hi Shri, Thanks for your message and all the helpful tips. If the TSVISetVariableBounds() functions are available now, I would like to try them as well. Is the interface essentially the same as with SNESVISetVariableBounds()? I will get back to you and Barry when I have had a chance to modify my application. For my problem, I believe it makes sense to have bounds on one of the variables as well as one function of the variables. The two relevant degrees of freedom are the block voltage (one for each finite difference block) and the electrode voltage (one for each electrode, which may be present in multiple blocks). The electrode voltage should keep a constant phase while the magnitude is constrained between zero and some maximum. The block voltages near the electrodes depend on the electrode voltages as well as the neighbouring block voltages. The electrode current for a given electrode is a function of its electrode voltage and several nearby block voltages, and should be constrained in magnitude between zero and some maximum (and the phase unconstrained). Would I need to use SNESVISetComputeVariableBounds since the electrode current is a function of the other variables? Would any other provisions need to be made for the block voltages, since they depend on the electrode voltages? ? You cannot directly make a bound on some function of other ? variables. Instead you need to introduce another variable that is ? defined to be equal to that function and put the bound on that new ? variable. ? Barry Thank you again, Jon On 2011-07-25, at 11:58 AM, Shri wrote: - Original Message - On Jul 22, 2011, at 4:16 PM, Jonathan Backs wrote: Barry, Thank you so much for your response. Lucky, indeed! I look forward to trying out these new features. I neglected to mention in my original post that my electrical problem is part of a DAE, which includes a time-dependent heating problem. Can SNESVI constraints be used in conjunction with TSSetIFunction() and TSSetIJacobian() as well? (Of course, I only need the constraints for the time-independent electrical portion.) We have not yet put that in but Shri is starting to look at that just now. Basically we would have a TSVISetVariableBounds() and handle everything from there. I suggest you start with a simple time electrical portion with constraints to explore and we'll go from there. Shri is actually a electrical networks guy and so can speak your language. ? I've added TSVISetVariableBounds() for setting the bounds on the ? variables using the TS object directly. A few things that i want to mention here about using the variational inequality nonlinear solver (SNESVI). i) Use the runtime option -snes_type vi or explicitly set SNESSetType(snes,SNESVI). ii) There are two tested algorithms currently available, (a) semismooth (-snes_vi_type ss) and (b) active set or reduced space (-snes_vi_type rs). iii) Take a look at example,ex61.c, in src/snes/examples/tutorials which is a time-stepping nonlinear problem with constraints on the variables. This example does not use the TS object directly but rather a time-loop is explicitly written. Another example,ex8.c, in src/snes/examples/tests/ is a minimum surface area problem which uses SNESVI. By the way, does your problem have bounds on the variables
[petsc-users] Conditional Constraints
- Original Message - On Jul 22, 2011, at 4:16 PM, Jonathan Backs wrote: Barry, Thank you so much for your response. Lucky, indeed! I look forward to trying out these new features. I neglected to mention in my original post that my electrical problem is part of a DAE, which includes a time-dependent heating problem. Can SNESVI constraints be used in conjunction with TSSetIFunction() and TSSetIJacobian() as well? (Of course, I only need the constraints for the time-independent electrical portion.) We have not yet put that in but Shri is starting to look at that just now. Basically we would have a TSVISetVariableBounds() and handle everything from there. I suggest you start with a simple time electrical portion with constraints to explore and we'll go from there. Shri is actually a electrical networks guy and so can speak your language. I've added TSVISetVariableBounds() for setting the bounds on the variables using the TS object directly. A few things that i want to mention here about using the variational inequality nonlinear solver (SNESVI). i) Use the runtime option -snes_type vi or explicitly set SNESSetType(snes,SNESVI). ii) There are two tested algorithms currently available, (a) semismooth (-snes_vi_type ss) and (b) active set or reduced space (-snes_vi_type rs). iii) Take a look at example,ex61.c, in src/snes/examples/tutorials which is a time-stepping nonlinear problem with constraints on the variables. This example does not use the TS object directly but rather a time-loop is explicitly written. Another example,ex8.c, in src/snes/examples/tests/ is a minimum surface area problem which uses SNESVI. By the way, does your problem have bounds on the variables or bounds on some function of the variables? Shri Barry Thank you again, Jon On 2011-07-22, at 2:46 PM, Barry Smith wrote: Jon, You may be in luck. In PETSc-dev http://www.mcs.anl.gov/petsc/petsc-as/developers/index.html we have now implemented variational inequality non-linear solvers with box constraints. That is one has a set of variables u and algebraic equations F(u) = 0 (say of size n each) but in addition one has constraints lu = u = uu where some or all of lu may be negative infinity (no constraints) and some or all of uu may be infinity (no constraints). There is also a constraint on the sign of F() for those equations associated with active constraints. If your constraints are not in this form sometimes you can introduce new additional variables to get it in this form. Read up a little on variational inequalities on the web. To use this you provide the usual SNES function and Jacobian but you also provide SNESVISetVariableBounds() there is a manual page for this function and for for SNESVI at http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-dev/docs/index.html (the link is broken right now but Satish is fixing it). Plus several examples in src/snes/examples/tutorials (in petsc-dev). This is all new code so we would be interested in your feedback. Barry On Jul 22, 2011, at 3:33 PM, Jonathan Backs wrote: Hi, I am trying to add a constraint feature to my first PETSc application, which uses the finite difference method to calculate the potential distribution produced by a collection of electrodes in a resistive medium. I would like to make this simulation more realistic by imposing a maximum electric current and a maximum potential difference that can be supplied to each electrode by the power supply. If the medium between the electrodes is very conductive, the current maximum would be exceeded by the maximum potential difference, so the potential difference should be decreased from maximum until it produces the maximum current. On the other hand, the potential difference between the electrodes should remain at maximum as long as the current remains below maximum (say, for a less conductive medium). I added an extra degree of freedom (the electrode voltages) to my DMDA, and I developed a set of conditional expressions that describe the above constraints, but one problem is that the logic relies on if-then-else decisions that are made when forming the function/residual and the Jacobian. Once these decisions are made, of course, the conditions are not checked again until the next function or Jacobian evaluation. The non-linear solver then tends to oscillate between extreme solutions to the opposing conditions with each iteration, and never converges towards a reasonable solution. Is there a better strategy for solving such problems? Does PETSc offer mechanisms to aid in their solution? I would very much appreciate any hints. Thank you for your time, Jon
[petsc-users] Error for ilu
Hi, I am getting this error for ilu. Did anything get changed for ilu? I get the same error for debug mode as well. ./TS -pc_type ilu [0]PETSC ERROR: - Error Message [0]PETSC ERROR: No support for this operation for this object type! [0]PETSC ERROR: Matrix type seqaij symbolic ILU! [0]PETSC ERROR: [0]PETSC ERROR: Petsc Development HG revision: e36155ef11881fd94f40c202d1aacce6d908f14b HG Date: Mon Apr 25 11:24:07 2011 -0500 [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: ./TS on a opt-mode named karlin.math.iit.edu by abhyshr Tue Jun 21 01:09:20 2011 [0]PETSC ERROR: Libraries linked from /home/abhyshr/petsc-dev/opt-mode/lib [0]PETSC ERROR: Configure run at Mon Apr 25 18:22:34 2011 [0]PETSC ERROR: Configure options --download-blacs --download-hypre --download-mpich --download-mumps --download-parmetis --download-plapack --download-scalapack --download-scotch --download-superlu --download-superlu_dist --download-umfpack --with-cc=gcc --with-cxx=g++ --with-debugging=0 --with-fc=gfortran --with-clanguage=cxx --download-f-blas-lapack [0]PETSC ERROR: [0]PETSC ERROR: MatILUFactorSymbolic() line 6011 in src/mat/interface/matrix.c [0]PETSC ERROR: PCSetUp_ILU() line 216 in src/ksp/pc/impls/factor/ilu/ilu.c [0]PETSC ERROR: PCSetUp() line 819 in src/ksp/pc/interface/precon.c [0]PETSC ERROR: KSPSetUp() line 261 in src/ksp/ksp/interface/itfunc.c [0]PETSC ERROR: KSPSolve() line 383 in src/ksp/ksp/interface/itfunc.c [0]PETSC ERROR: SNES_KSPSolve() line 3124 in src/snes/interface/snes.c [0]PETSC ERROR: SNESSolve_LS() line 190 in src/snes/impls/ls/ls.c [0]PETSC ERROR: SNESSolve() line 2424 in src/snes/interface/snes.c [0]PETSC ERROR: main() line 160 in TS-dir/TS.c application called MPI_Abort(MPI_COMM_WORLD, 56) - process 0 [unset]: aborting job: application called MPI_Abort(MPI_COMM_WORLD, 56) - process 0 Thanks, Shri
[petsc-users] Error for ilu
Sorry, my mistake..since 3 other people are working on this project someone added pc_factor_package mumps in the options file and hence ilu was giving that error. Thanks, Shri - Original Message - Shri, I notice this is on Karlin. Does the same code work on other machines? Does petsc examples work on this machine? Hong On Tue, Jun 21, 2011 at 4:12 AM, Jed Brown jed at 59a2.org wrote: On Tue, Jun 21, 2011 at 08:15, Shri abhyshr at mcs.anl.gov wrote: I am getting this error for ilu. Did anything get changed for ilu? I get the same error for debug mode as well. ./TS -pc_type ilu Can you reproduce with an example? I haven't noticed any problems. Maybe you have memory corruption or something else unexpected causing MatGetFactor to not be called. Maybe try setting a breakpoint at ilu.c:213.
[petsc-users] Error for ilu
The error message is indeed confusing. Perhaps a better error message would be 'Cannot use -pc_type ilu with mumps'. - Original Message - Hmm, then the error message is not complete enough. It printed seqaij but nothing about mumps hence confused everyone who read your email message. So how come mumps doesn't come up somewhere in the error message? Barry On Jun 21, 2011, at 11:14 AM, Shri wrote: Sorry, my mistake..since 3 other people are working on this project someone added pc_factor_package mumps in the options file and hence ilu was giving that error. Thanks, Shri - Original Message - Shri, I notice this is on Karlin. Does the same code work on other machines? Does petsc examples work on this machine? Hong On Tue, Jun 21, 2011 at 4:12 AM, Jed Brown jed at 59a2.org wrote: On Tue, Jun 21, 2011 at 08:15, Shri abhyshr at mcs.anl.gov wrote: I am getting this error for ilu. Did anything get changed for ilu? I get the same error for debug mode as well. ./TS -pc_type ilu Can you reproduce with an example? I haven't noticed any problems. Maybe you have memory corruption or something else unexpected causing MatGetFactor to not be called. Maybe try setting a breakpoint at ilu.c:213.
[petsc-users] Error for ilu
Thanks for the patch. I tried to run the code but it spews an error Undefined symbols: MatMFFDSetBase_MFFD(_p_Mat*, _p_Vec*, _p_Vec*), referenced from: MatAssemblyEnd_SNESMF(_p_Mat*, MatAssemblyType) in libpetsc.a(snesmfj.o) _MatMFFDSetBase_SNESMF in libpetsc.a(snesmfj.o) make test is also giving the same error - Original Message - On Tue, Jun 21, 2011 at 20:02, Shri abhyshr at mcs.anl.gov wrote: The error message is indeed confusing. Perhaps a better error message would be 'Cannot use -pc_type ilu with mumps'. Did you try with the patch I pushed 45 minutes ago? [0]PETSC ERROR: No support for this operation for this object type! [0]PETSC ERROR: Matrix type seqaij symbolic ILU using solver package mumps! -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110621/563e51fc/attachment.htm
[petsc-users] How to get a duoble array into petsc Vec?
You could use VecCreateSeqWithArray() http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Vec/VecCreateSeqWithArray.html Shri - Original Message - Hi, I am trying to put data form a double c array to an petsc vector. I tried to loop over each element like this petsc Vec[i] = double array[i], but I need a cast for that?.has anybody an idea? Thanks J?rgen -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110519/860007f8/attachment.htm
[petsc-users] questions about parallel vectors
- Original Message - Hello everyone, I am trying to create a sparse matrix A, which depends on a parallel vector v. For example, my function looks like this: Mat myfun(MPI_Comm comm, Vec v, other parameters). When I set the value A(i,j) = v[k], v[k] may not be obtained by VecGetValues since that operation can only get values on the same processor. I am thinking 1) create v as an array and pass this array into myfun. 2) create another vector v2, which is a full copy of parallel v through VecScatter. Do you need all the vector elements on each processor to set the matrix values or just a subset of them? ' If you need all the vector elements then you can use VecScattertoAll http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Vec/VecScatterCreateToAll.html If you only need a subset then you could create v as a ghosted vector. See the example http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/src/vec/vec/examples/tutorials/ex9.c.html 3) when I first create the initial vec v, using VecCreate(PETSC_COMM_SELF,v) or VecCreateSeq. Does this guarantee that all the processors creating matrix A have all the components of vector v? I think 1) and 2) are going to work, but not sure about option 3). I have no idea which would have better performance. Can you give me some suggestions on how to handle this problem? Thanks. Another quick question, what is the difference between PetscViewerSetFormat and PetscViewerPushFormat? Best, Xiangdong
[petsc-users] Least squares: using unpreconditioned LSQR
Use the option -pc_type none. - Original Message - Hi, I wanted to solve some least squares problems using PETSc. My test matrix is size 3x2 but I wish to use this code for solving large ill-conditioned rectangular systems later. Looking at the PETSc manual I the found the KSPLSQR routine which implements the LSQR algorithm. However I am unsure how to use this routine. I am pasting the lines of the code which I use to set up the solver. Through the terminal I pass the option -ksp_type lsqr while running exectuable. ierr = KSPCreate(PETSC_COMM_WORLD,ksp);CHKERRQ(ierr); ierr = KSPSetOperators(ksp,A,PETSC_NULL,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); As you can see I have used PETSC_NULL for the preconditioner matrix since I wish to use the *unpreconditioned* version of the LSQR algorithm. This gives me an error. If I pass the matrix A it gives me an error again. I am not sure how to tell PETSc not to use a preconditioner. Could you pease tell me how I should use KSPSetOperators statement in this case to use the unpreconditioned algorithm. If you have a better sparse matrix least squares algorithm implemented please let me know. -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110217/cc7a6d20/attachment.htm
[petsc-users] Trouble with solving 5x5 system... Ay=x on 2 or more processors
You are not setting the 'global' and the 'local' sizes of vector x correctly as the error message says. The global size of the matrix is 5X5 while the global size of the vector x is 10 !!! - Original Message - Hi, I am trying to solve a linear system where Ay=x and A is a 5x5 matrix stored in a binary file called 'square' and x=[1;1;1;1;1] I am trying to display the matrix A, and the vectors x(rhs) and y(soln) in that order to standard output. On running my code on a single processor the answer returned is accurate. But on using 2 processors I get weird error messages PART of which says [0]PETSC ERROR: [1]PETSC ERROR: - Error Message [1]PETSC ERROR: Nonconforming object sizes! [1]PETSC ERROR: Mat mat,Vec x: global dim 5 10! Also somehow the vector x gets displayed TWICE when run on two processes. A however gets displayed ONCE (as it should!!) I am attaching the output I get when I run on 1 process and the when I run the same code on 2 processes. please let me know where I could be going wrong. (1) This is what I get on running on ONE process: (Here system is solved successfully) gaurish108 at gaurish108-laptop:~/Desktop$ $PETSC_DIR/$PETSC_ARCH/bin/mpiexec -n 1 ./ex4 -f square 1.5761308167754828e-01 1.4188633862721534e-01 6.5574069915658684e-01 7.5774013057833345e-01 7.0604608801960878e-01 9.7059278176061570e-01 4.2176128262627499e-01 3.5711678574189554e-02 7.4313246812491618e-01 3.1832846377420676e-02 9.5716694824294557e-01 9.1573552518906709e-01 8.4912930586877711e-01 3.9222701953416816e-01 2.7692298496088996e-01 4.8537564872284122e-01 7.9220732955955442e-01 9.3399324775755055e-01 6.5547789017755664e-01 4.6171390631153941e-02 8.002804680011e-01 9.5949242639290300e-01 6.7873515485777347e-01 1.7118668781156177e-01 9.7131781235847536e-02 Process [0] 1 1 1 1 1 KSPGetIterationNumber 5 KSPGetResidualNorm 0.00 Process [0] -0.810214 2.33178 -1.31131 1.09323 1.17322 gaurish108 at gaurish108-laptop:~/Desktop$ % This is what I get on running on TWO processes: (2) gaurish108 at gaurish108-laptop:~/Desktop$ $PETSC_DIR/$PETSC_ARCH/bin/mpiexec -n 2 ./ex4 -f square 1.5761308167754828e-01 1.4188633862721534e-01 6.5574069915658684e-01 7.5774013057833345e-01 7.0604608801960878e-01 9.7059278176061570e-01 4.2176128262627499e-01 3.5711678574189554e-02 7.4313246812491618e-01 3.1832846377420676e-02 9.5716694824294557e-01 9.1573552518906709e-01 8.4912930586877711e-01 3.9222701953416816e-01 2.7692298496088996e-01 4.8537564872284122e-01 7.9220732955955442e-01 9.3399324775755055e-01 6.5547789017755664e-01 4.6171390631153941e-02 8.002804680011e-01 9.5949242639290300e-01 6.7873515485777347e-01 1.7118668781156177e-01 9.7131781235847536e-02 Process [0] 1 1 1 1 1 Process [1] 1 1 1 1 1 [0]PETSC ERROR: [1]PETSC ERROR: - Error Message [1]PETSC ERROR: Nonconforming object sizes! [1]PETSC ERROR: Mat mat,Vec x: global dim 5 10! [1]PETSC ERROR: [1]PETSC ERROR: Petsc Release Version 3.1.0, Patch 5, Mon Sep 27 11:51:54 CDT 2010 [1]PETSC ERROR: See docs/changes/index.html for recent updates. [1]PETSC ERROR: See docs/faq.html for hints about trouble shooting. [1]PETSC ERROR: - Error Message [0]PETSC ERROR: Nonconforming object sizes! [0]PETSC ERROR: Mat mat,Vec x: global dim 5 10! [0]PETSC ERROR: See docs/index.html for manual pages. [1]PETSC ERROR: [1]PETSC ERROR: ./ex4 on a linux-gnu named gaurish108-laptop by gaurish108 Mon Jan 17 23:49:18 2011 [1]PETSC ERROR: Libraries linked from /home/gaurish108/Desktop/ResearchMeetings/SUPERPETS/petsc-3.1-p5/linux-gnu-c-debug/lib [1]PETSC ERROR: Configure run at Sat Nov 13 20:34:38 2010 [1]PETSC ERROR: Configure options --with-cc=gcc --with-fc=gfortran --download-f-blas-lapack=1 --download-mpich=1 --download-superlu_dist=1 --download-parmetis=1 --with-superlu_dist=1 --with-parmetis=1 [1]PETSC ERROR: [1]PETSC ERROR: MatMultTranspose() line 1947 in src/mat/interface/matrix.c [1]PETSC ERROR: KSPSolve_CGNE() line 103 in src/ksp/ksp/impls/cg/cgne/cgne.c [1]PETSC ERROR: KSPSolve() line 396 in src/ksp/ksp/interface/itfunc.c [1]PETSC ERROR: main() line 78 in src/mat/examples/tutorials/ex4.c application called MPI_Abort(MPI_COMM_WORLD, 60) - process 1[cli_1]: aborting job: application called MPI_Abort(MPI_COMM_WORLD, 60) - process 1 [0]0:Return code = 0, signaled with Interrupt [0]1:Return code = 60 gaurish108 at gaurish108-laptop:~/Desktop$ -- next part -- An HTML attachment
[petsc-users] Reading matrices into PETSc
Find attached a routine which reads matrix data from an ASCII file in i j value format and creates a seqaij matrix. - Original Message - I am sorry, but I am not yet clear on how to do this. I read ex32.c and ex72.c but I am still confused. What is ASCII 'slap' format? How should matrix be supplied to PETSc? My matrix is a 2000x1900 matrix given in the format MATLAB stores sparse matrices. i.e [row, column, non-zero-entry] format. On Tue, Dec 21, 2010 at 11:03 PM, Gaurish Telang gaurish108 at gmail.com wrote: i have this large text file containing a matrix. This text file contains the non-zero entries of a very large sparse matrix The first two columns indicate the position of the non-zero entry and the last column the actual non-zero value it self for example the matrix 1 0 8 0 0 5 6 0 0 is written in the text file in the form of 1 1 1 1 3 8 2 3 5 3 1 6 This is the standard [(row ,column), non-zero] entry format. i want PETSc to load this matrix from the text file i am not sure how to do that. What commands do I use? I am new to PETSc, so some detail in the explanation will be really helpful. Sincere thanks, Gaurish. -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20101222/0ddb21ee/attachment-0001.htm -- next part -- A non-text attachment was scrubbed... Name: ReadMatFromFile.c Type: text/x-csrc Size: 2125 bytes Desc: not available URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20101222/0ddb21ee/attachment-0001.c
[petsc-users] Bug when multipling a SBAIJ matrix with a vector
Fixed and pushed to petsc-3.1 - Barry Smith bsmith at mcs.anl.gov wrote: I've asked Hong and Shri to fix the bug I introduced that sometimes makes it require a diagonal. Once the bug is fixed 3.1 will allow empty diagonals again. Barry On Oct 21, 2010, at 8:08 AM, Jed Brown wrote: On Thu, Oct 21, 2010 at 14:41, Andreas Hauffe andreas.hauffe at tu-dresden.de wrote: I get no error and this example delivers the right result for PETSC 3.0. What did change from 3.0 to 3.1? Many of the matrix kernels were optimized, this one was changed in relatively simple ways, but the data structures for factorization were changed completely. In this case, Barry's name is on the commit message http://petsc.cs.iit.edu/petsc/petsc-dev/rev/0308bd570415#l2.33 He should have put that behavioral change in the 3.1 changelog. Or maybe the cost of handling empty rows is actually not significant so it should still be supported. I can't tell from the commit message, but empty rows are still supported by the code in MatMult_SeqSBAIJ_2, so it probably should work for block size 1 as well. Do you really not get an error message when you run your example code with 3.1 built with debug support? I just ran your code with 3.1 and I get a (bad) error message. Jed
[petsc-users] matrix-matrix addition
Ok. I'll work on it this afternoon. - Barry Smith bsmith at mcs.anl.gov wrote: On Sep 27, 2010, at 12:04 PM, Jed Brown wrote: On Mon, Sep 27, 2010 at 19:01, Barry Smith bsmith at mcs.anl.gov wrote: If you do not merge the column indices from A and B but simply count them all you will get over allocation, at worst a factor of two, though usually it would just be a little. You count the number of *unique* entries. You have two indices that run through the arrays aj and bj, counting duplicates only once. This is easy to do since they are both sorted. Yes, sounds like it could work and definitely better than my suggestion. Shri, can you try this? Barry Jed
[petsc-users] Fwd: Re: Question on zero sized DAs
Barry, Sorry i forgot to mention that i'm using petsc-3.1 and not petsc-dev. I'll add your change to petsc-3.1. Thanks, Shri - Barry Smith bsmith at mcs.anl.gov wrote: Pull and recompile in src/dm/da/src/ then try running again and report any problems to petsc-maint Do you have a stencil width greater than 0? If you are using a stencil width of zero we can try to bypass this error check and see if things go through. BTW: this is a petsc-maint question :-) or petsc-users Barry On Sep 8, 2010, at 8:56 PM, Shri wrote: I have three DAs da1,da2, and da3 (1-dimensional DAs) for my power system project where da1 is for the generator subsystem, da2 is for the network subsystem, and da3 is for the load subsystem. Here,the load subsystem size is zero and hence i need to create da3 with zero nodes. However,DACreate1D does not allow zero nodes or zero degrees of freedom and gives an error [0]PETSC ERROR: - Error Message [0]PETSC ERROR: Argument out of range! [0]PETSC ERROR: More processors than data points! 1 0! [0]PETSC ERROR: [0]PETSC ERROR: Petsc Release Version 3.1.0, Patch 4, unknown [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: /Users/Shri/TS-dev/TS1 on a osx-debug named dhcp55.swwn2.iit.edu by Shri Wed Sep 8 19:56:14 2010 [0]PETSC ERROR: Libraries linked from /Users/Shri/petsc-3.1/osx-debug/lib [0]PETSC ERROR: Configure run at Fri Jul 23 16:06:56 2010 [0]PETSC ERROR: Configure options --with-mpi-dir=/Users/Shri/packages/install/mpich2 --with-clanguage=cxx [0]PETSC ERROR: [0]PETSC ERROR: DACreate_1D() line 152 in src/dm/da/src/da1.c [0]PETSC ERROR: DASetType() line 48 in src/dm/da/src/dareg.c [0]PETSC ERROR: DASetTypeFromOptions_Private() line 65 in src/dm/da/src/dacreate.c [0]PETSC ERROR: DASetFromOptions() line 131 in src/dm/da/src/dacreate.c [0]PETSC ERROR: DACreate1d() line 393 in src/dm/da/src/da1.c [0]PETSC ERROR: CreateLoadSysDA() line 80 in createsubsysda.c [0]PETSC ERROR: main() line 64 in TS1.c application called MPI_Abort(MPI_COMM_WORLD, 63) - process 0 1) Is it possible to create 1d DAs with zero nodes?? 2) I can try using vectors instead of DAs for the three subsystems since zero sized vectors can be created. However, i want to pack these three subsystem vectors later on similar to what the DMComposite object does. Did PETSc have a vector packer object similar to DMComposite in any of the previous releases? Thanks, Shri
