Thank you Lisandro,<br><br>I like the idea of constructing Schur complement of Schur complement matrix. I will give it a try.<br><br>Thanks again for the suggestions and the link <br><br>Waad :)<br><br><b><i>Lisandro Dalcin &lt;dalcinl@gmail.com&gt;</i></b> wrote:<blockquote class="replbq" style="border-left: 2px solid rgb(16, 16, 255); margin-left: 5px; padding-left: 5px;"> On 5/23/08, Barry Smith <bsmith@mcs.anl.gov> wrote:<br>&gt; I notice Lapack solver for the interior<br>&gt; problem takes a lot of time compare to the iterative solver for the<br>&gt; interface, so now I am replacing the direct factorization with petsc KSP<br>&gt; solver.<br>&gt; &gt;<br>&gt;<br>&gt;    In my experience you will want to use a KSP direct solve, for example<br>&gt; just -ksp_type preonly -pc_type lu; using an iterative solver in such Schur<br>&gt; complement is<br>&gt;  always way to slow.<br><br>Indeed!!<br><br>Waad, do not do that. If your interior problem is big, the only way to<br>go
 is to implement subdomain subpartitioning. That is, in your local<br>subdomain, you have to select some DOF's and 'mark' them as<br>'interface' nodes. Then you end-up having many A_ii's at a local<br>subdomain, each corresponding with a sub-subdomain... Of course,<br>implementing all this logic is not trivial at all..<br><br>My implementation of all this stuff can be found in petsc4py SVN<br>repository hosted at Google Code. Here you have the (long) link:<br><br>http://petsc4py.googlecode.com/svn/trunk/petsc/lib/ext/petsc/src/ksp/pc/impls/<br><br><br>In directory 'schur' you have a version requiring MATIS matrix type<br>(this is somewhat natural in the context of substructuring and finite<br>elements methods). This corresponds to what Y.Sadd's book calls<br>'edge-based' partitionings.<br><br>In directory 'schur_aij' hou have a version (never carefully tested)<br>working with MATMPIAIJ matrices. This corresponds to what Y.Sadd's<br>book calls 'vertex-based' partitionings
 (more typical to appear in<br>finite difference methods).<br><br><br><br><br><br><br>&gt; &gt; I would like very much to have a look at your implementation, and I think<br>&gt; that will be very useful to me.<br>&gt; &gt;<br>&gt; &gt; Thanks<br>&gt; &gt;<br>&gt; &gt; Waad<br>&gt; &gt;<br>&gt; &gt; Lisandro Dalcin <dalcinl@gmail.com> wrote: On 5/20/08, Waad Subber wrote:<br>&gt; &gt; &gt; 1) How do you actually get the local Schur complements. You<br>&gt; &gt; &gt; explicitelly compute its entries, or do you compute it after computing<br>&gt; &gt; &gt; the inverse (or LU factors) of a 'local' matrix?<br>&gt; &gt; &gt;<br>&gt; &gt; &gt; I construct the local Schur complement matrices after getting the<br>&gt; inversion<br>&gt; &gt; &gt; of A_II matrix for each subdomain.<br>&gt; &gt;<br>&gt; &gt; Fine,<br>&gt; &gt;<br>&gt; &gt; &gt; 2) Your R_i matrix is actually a matrix? In that case, it is a trivial<br>&gt; &gt; &gt; restrinction operation with ones and zeros? Or R_i is
 actually a<br>&gt; &gt; &gt; VecScatter?<br>&gt; &gt; &gt;<br>&gt; &gt; &gt; R_i is the restriction matrix maps the global boundary nodes to the<br>&gt; local<br>&gt; &gt; &gt; boundary nodes and its entries is zero and one I store it as spare<br>&gt; matrix,<br>&gt; &gt; &gt; so only I need to store the nonzero entries which one entry per a row<br>&gt; &gt;<br>&gt; &gt; I believe a VecScatter will perform much better for this task.<br>&gt; &gt;<br>&gt; &gt;<br>&gt; &gt; &gt; And finally: are you trying to apply a Krylov method over the global<br>&gt; &gt; &gt; Schur complement? In such a case, are you going to implement a<br>&gt; &gt; &gt; preconditioner for it?<br>&gt; &gt; &gt;<br>&gt; &gt; &gt; Yes, that what I am trying to do<br>&gt; &gt;<br>&gt; &gt; Well, please let me make some comments. I've spent many days and month<br>&gt; &gt; optimizing Schur complement iterations, and I ended giving up. I was<br>&gt; &gt; never able to get it perform better than ASM
 preconditioner (iff<br>&gt; &gt; appropriatelly used, ie. solving local problems with LU, and<br>&gt; &gt; implementing subdomain subpartitioning the smart way, not the way<br>&gt; &gt; currently implemented in PETSc, were subpartitioning is done by chunks<br>&gt; &gt; of continuous rows).<br>&gt; &gt;<br>&gt; &gt; If you are doing research on this, I would love to know your<br>&gt; &gt; conclusion when you get your work done. If you are doing all this just<br>&gt; &gt; with the hope of getting better running times, well, remember my above<br>&gt; &gt; comments but also remember that I do not consider myself a smart guy<br>&gt; &gt; ;-)<br>&gt; &gt;<br>&gt; &gt; As I said before, I worked hard for implementing general Schur<br>&gt; &gt; complement iteration. All this code is avalable in the SVN repository<br>&gt; &gt; of petsc4py (PETSc for Python), but it could be easily stripped out<br>&gt; &gt; for use in any PETSc-based code in C/C++. This implementation
 requires<br>&gt; &gt; the use of a MATIS matrix type (there is also a separate<br>&gt; &gt; implementation for MATMPIAIJ maatrices), I've implemented subdomain<br>&gt; &gt; subpartitioning (using a simple recursive graph splitting procedure<br>&gt; &gt; reusing matrix reordering routines built-in in PETSc, could be done<br>&gt; &gt; better with METIS); when the A_ii problems are large, their LU<br>&gt; &gt; factorization can be a real bootleneck. I've even implemented a<br>&gt; &gt; global preconditioner operation for the interface problem, based on<br>&gt; &gt; iterating over a 'strip' of nodes around the interface; it improves<br>&gt; &gt; convergence and is usefull for ill-conditioned systems, but the costs<br>&gt; &gt; are increased.<br>&gt; &gt;<br>&gt; &gt; If you ever want to take a look at my implemention for try to use it,<br>&gt; &gt; or perhaps take ideas for your own implementation, let me know.<br>&gt; &gt;<br>&gt; &gt;<br>&gt; &gt;<br>&gt; &gt;<br>&gt;
 &gt;<br>&gt; &gt; &gt; &gt; Now having the Schur complement matrix for each subdomain, I need to<br>&gt; solve<br>&gt; &gt; &gt; &gt; the interface problem<br>&gt; &gt; &gt; (Sum[R_i^T*S_i*R_i])u=Sum[R_i^T*g_i],<br>&gt; &gt; &gt; &gt; .. i=1.. to No. of process (subdomains) in parallel.<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; For the global vector I construct one MPI vector and use VecGetArray<br>&gt; ()<br>&gt; &gt; &gt; for<br>&gt; &gt; &gt; &gt; each of the sequential vector then use VecSetValues () to add the<br>&gt; values<br>&gt; &gt; &gt; &gt; into the global MPI vector. That works fine.<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; However for the global schur complement matix I try the same idea by<br>&gt; &gt; &gt; &gt; creating one parallel MPIAIJ matrix and using MatGetArray( ) and<br>&gt; &gt; &gt; &gt; MatSetValues () in order to add the values to the global matrix.<br>&gt; &gt; &gt; &gt; MatGetArray( ) gives me only the values without indices, so I
 don't<br>&gt; know<br>&gt; &gt; &gt; how<br>&gt; &gt; &gt; &gt; to add these valuse to the global MPI matrix.<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; Thanks agin<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; Waad<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; Barry Smith wrote:<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; On May 20, 2008, at 3:16 PM, Waad Subber wrote:<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; Thank you Matt,<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; Any suggestion to solve the problem I am trying to tackle. I want to<br>&gt; &gt; &gt; &gt; &gt; solve a linear system:<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; Sum(A_i) u= Sum(f_i) , i=1.... to No. of CPUs.<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; Where A_i is a sparse sequential matrix and f_i is a sequential<br>&gt; &gt; &gt; &gt; &gt; vector. Each CPU has one matrix and one vector of the same size. Now<br>&gt; &gt; &gt; &gt; &gt; I want to sum up
 and solve the system in parallel.<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; Does each A_i have nonzero entries (mostly) associated with one<br>&gt; &gt; &gt; &gt; part of the matrix? Or does each process have values<br>&gt; &gt; &gt; &gt; scattered all around the matrix?<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; In the former case you should simply create one parallel MPIAIJ<br>&gt; &gt; &gt; &gt; matrix and call MatSetValues() to put the values<br>&gt; &gt; &gt; &gt; into it. We don't have any kind of support for the later case, perhaps<br>&gt; &gt; &gt; &gt; if you describe how the matrix entries come about someone<br>&gt; &gt; &gt; &gt; would have suggestions on how to proceed.<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; Barry<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; Thanks again<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; Waad<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt;
 Matthew Knepley wrote: On Tue, May 20, 2008 at<br>&gt; &gt; &gt; &gt; &gt; 2:12 PM, Waad Subber wrote:<br>&gt; &gt; &gt; &gt; &gt; &gt; Hi,<br>&gt; &gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; &gt; I am trying to construct a sparse parallel matrix (MPIAIJ) by<br>&gt; &gt; &gt; &gt; &gt; adding up<br>&gt; &gt; &gt; &gt; &gt; &gt; sparse sequential matrices (SeqAIJ) from each CPU. I am using<br>&gt; &gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; &gt; MatMerge_SeqsToMPI(MPI_Comm comm,Mat seqmat,PetscInt m,PetscInt<br>&gt; &gt; &gt; &gt; &gt; n,MatReuse<br>&gt; &gt; &gt; &gt; &gt; &gt; scall,Mat *mpimat)<br>&gt; &gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; &gt; to do that. However, when I compile the code I get the following<br>&gt; &gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; &gt; undefined reference to `matmerge_seqstompi_'<br>&gt; &gt; &gt; &gt; &gt; &gt; collect2: ld returned 1 exit status<br>&gt; &gt; &gt; &gt; &gt; &gt; make: *** [all] Error
 1<br>&gt; &gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; &gt; Am I using this function correctly ?<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; These have no Fortran bindings right now.<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; Matt<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; &gt; Thanks<br>&gt; &gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; &gt; Waad<br>&gt; &gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt; --<br>&gt; &gt; &gt; &gt; &gt; What most experimenters take for granted before they begin their<br>&gt; &gt; &gt; &gt; &gt; experiments is infinitely more interesting than any results to which<br>&gt; &gt; &gt; &gt; &gt; their experiments lead.<br>&gt; &gt; &gt; &gt; &gt; -- Norbert Wiener<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt;
 &gt;<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt; &gt;<br>&gt; &gt; &gt;<br>&gt; &gt; &gt;<br>&gt; &gt; &gt; --<br>&gt; &gt; &gt; Lisandro Dalcín<br>&gt; &gt; &gt; ---------------<br>&gt; &gt; &gt; Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)<br>&gt; &gt; &gt; Instituto de Desarrollo Tecnológico para la Industria Química (INTEC)<br>&gt; &gt; &gt; Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)<br>&gt; &gt; &gt; PTLC - Güemes 3450, (3000) Santa Fe, Argentina<br>&gt; &gt; &gt; Tel/Fax: +54-(0)342-451.1594<br>&gt; &gt; &gt;<br>&gt; &gt; &gt;<br>&gt; &gt; &gt;<br>&gt; &gt; &gt;<br>&gt; &gt; &gt;<br>&gt; &gt;<br>&gt; &gt;<br>&gt; &gt; --<br>&gt; &gt; Lisandro Dalcín<br>&gt; &gt; ---------------<br>&gt; &gt; Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)<br>&gt; &gt; Instituto de Desarrollo Tecnológico para la Industria Química (INTEC)<br>&gt; &gt; Consejo Nacional de Investigaciones Científicas y Técnicas
 (CONICET)<br>&gt; &gt; PTLC - Güemes 3450, (3000) Santa Fe, Argentina<br>&gt; &gt; Tel/Fax: +54-(0)342-451.1594<br>&gt; &gt;<br>&gt; &gt;<br>&gt; &gt;<br>&gt; &gt;<br>&gt;<br>&gt;<br><br><br>-- <br>Lisandro Dalcín<br>---------------<br>Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)<br>Instituto de Desarrollo Tecnológico para la Industria Química (INTEC)<br>Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)<br>PTLC - Güemes 3450, (3000) Santa Fe, Argentina<br>Tel/Fax: +54-(0)342-451.1594<br><br></dalcinl@gmail.com></bsmith@mcs.anl.gov></blockquote><br><p>&#32;