Understood. Where are you putting the "few extra unknowns" in the vector and
matrix? On the first process, on the last process, some places in the middle of
the matrix?
We don't have any trivial code for copying a big matrix into a even larger
matrix directly because we frown on doing that. It is very wasteful in time and
memory.
The simplest way to do it is call MatGetRow() twice for each row, once to
get the nonzero locations for each row to determine the numbers needed for
preallocation and then the second time after the big matrix has been
preallocated to get the nonzero locations and numerical values for the row to
call MatSetValues() with to set that row into the bigger matrix. Note of course
when you call MatSetValues() you will need to shift the rows and column
locations to take into account the new rows and columns in the bigger matrix.
If you put the "extra unknowns" at the every end of the rows/columns on the
last process you won't have to shift.
Note that B being dense really messes up chances for load balancing since
its rows are dense and take a great deal of space so whatever process gets
those rows needs to have much less of the mesh.
The correct long term approach is to have libmesh provide the needed
functionality (for continuation) for the slightly larger matrix directly so
huge matrices do not need to be copied.
I noticed that libmesh has some functionality related to continuation. I do
not know if they handle it by creating the larger matrix and vector and filling
that up directly for finite elements. If they do then you should definitely
take a look at that and see if it can be extended for your case (ignore the
continuation algorithm they may be using, that is not relevant, the question is
if they generate the larger matrices and if you can leverage this).
The ultimate hack would be to (for example) assign the extra variables to the
end of the last process and hack lib mesh a little bit so the matrix it creates
(before it puts in the numerical values) has the extra rows and columns, that
libmesh will not put the values into but you will. Thus you get libmesh to fill
up the true final matrix for its finite element problem (not realizing the
matrix is a little bigger then it needs) directly, no copies of the data
needed. But this is bit tricky, you'll need to combine libmesh's preallocation
information with yours for the final columns and rows before you have lib mesh
put the numerical values in. Double check if they have any support for this
first.
Barry
> On May 29, 2019, at 6:29 PM, Manav Bhatia <bhatiama...@gmail.com> wrote:
>
> Thanks, Barry.
>
> I am working on a FE application (involving bifurcation behavior) with
> libMesh where I need to solve the system of equations along with a few extra
> unknowns that are not directly related to the FE mesh. I am able to assemble
> the n x 1 residual (R_fe) and n x n Jacobian (J_fe ) from my code and
> libMesh provides me with the sparsity pattern for this.
>
> Next, the system of equations that I need to solve is:
>
> [ J_fe A ] { dX } = { R_fe }
> [ B C ] { dV } = {R_ext }
>
> Where, C is a dense matrix of size m x m ( m << n ), A is n x m, B is m x n,
> R_ext is m x 1. A, B and C are dense matrixes. This comes from the bordered
> system for my path continuation solver.
>
> I have implemented a solver using Schur factorization ( this is outside of
> PETSc and does not use the FieldSplit construct ). This works well for most
> cases, except when J_fe is close to singular.
>
> I am now attempting to create a monolithic matrix that solves the complete
> system.
>
> Currently, the approach I am considering is to compute J_fe using my libMesh
> application, so that I don’t have to change that. I am defining a new matrix
> with the extra non-zero locations for A, B, C.
>
> With J_fe computed, I am looking to copy its non-zero entries to this new
> matrix. This is where I am stumbling since I don’t know how best to get the
> non-zero locations in J_fe. Maybe there is a better approach to copy from
> J_fe to the new matrix?
>
> I have looked through the nested matrix construct, but have not given this a
> serious consideration. Maybe I should? Note that I don’t want to solve J_fe
> and C separately (not as separate systems), so the field-split approach will
> not be suitable here.
>
> Also, I am currently using MUMPS for all my parallel solves.
>
> I would appreciate any advice.
>
> Regards,
> Manav
>
>
>> On May 29, 2019, at 6:07 PM, Smith, Barry F. <bsm...@mcs.anl.gov> wrote:
>>
>>
>> Manav,
>>
>> For parallel sparse matrices using the standard PETSc formats the matrix is
>> stored in two parts on each process (see the details in MatCreateAIJ()) thus
>> there is no inexpensive way to access directly the IJ locations as a single
>> local matrix. What are you hoping to use the information for? Perhaps we
>> have other suggestions on how to achieve the goal.
>>
>> Barry
>>
>>
>>> On May 29, 2019, at 2:27 PM, Manav Bhatia via petsc-users
>>> <petsc-users@mcs.anl.gov> wrote:
>>>
>>> Hi,
>>>
>>> Once a MPI-AIJ matrix has been assembled, is there a method to get the
>>> nonzero I-J locations? I see one for sequential matrices here:
>>> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatGetRowIJ.html
>>> , but not for parallel matrices.
>>>
>>> Regards,
>>> Manav
>>>
>>>
>>
>
