Newbie question about synchronization

[Knut Erik Teigen]

On Thu, 2007-03-29 at 03:42 -0400, Diego Nehab wrote:
> Hi,
> 
> I am very new to Petsc, and I am not very knowledgeable in
> numerical methods. Therefore, please forgive my ignorance.
> 
> I need to solve sparse (at most 7 non-zero elements per row)
> overconstrained linear systems, in the least squares sense.
> The systems are very large (typically 3M x 500k), but have a
> reasonably strong diagonal component that apparently makes
> them solvable.
> 
> In the past I have directly used the C version of
> C. C. Paige and M. A. Saunders LSQR library, with no problems.
> My interest in Petsc is the eternal quest for speed. That,
> I hope, could come in the form of better methods, a better matrix
> representation and, naturally, parallelism.
> 
> 1) Am I right to assume the only solver that accepts
> overconstrained systems is the lsqr?
> 
> So far, using only one process, everything is simple and
> works (it took me longer to compile and test MPI and Petsc
> then to write code that solves my problem :)).
> 
> When I move to multiple processes, the solver still works, but I
> couldn't figure out how to select one of the processes to
> consolidate the solution vector and use it to save a file on disk.
> I always get an error of the form
> 
>     [0]PETSC ERROR: Argument out of range!
>     [0]PETSC ERROR: Can only get local values, trying xxx!
> 
> I assume I must bring the solution vector back from the other
> processes, right?
> 
> 2) If so, how do I do this?
If you only want to save the results to disk, you can use the VecView
function. Just create a viewer, like e.g.
PetscViewerAsciiOpen(PETSC_COMM_WORLD,"filename",&viewer)
VecView(solution,viewer)
You can also output in binary format using BinaryOpen instead. Check the
chapter on Viewers in the manual.
If you need to gather the results on one processor for further
computations, I use standard MPI calls, like this(in Fortran):

 call VecGetArray(sol,sol_ptr,sol_i,ierr)

 ! gather solution on process 0
 if (rank==0) then 
   do i=low,high !copy local solution to global solution
     global_sol(i)=sol_ptr(sol_i+i)
   enddo
   do p=1,size-1 !recieve local solution from other processes
     call
MPI_Recv(local_sol,loc_n,MPI_REAL,p,1,PETSC_COMM_WORLD,istat,mpierr)
     do i=1,loc_n !copy local part to correct position in global
       global_sol(i+high)=local_sol(i)
       high=high+1
     enddo
   enddo
 else
  do j=1,loc_n
    local_sol(j)=sol_ptr(sol_i+j)
  enddo  
  call MPI_Send(local_sol,loc_n,MPI_REAL,0,1,PETSC_COMM_WORLD,mpierr)
 endif

 call VecRestoreArray(sol,sol_ptr,sol_i,ierr)

 !copy global solution vector back to grid array
 if (rank==0) then
   do j=1,jmax
     do i=1,imax
       T(i,j)=global_sol((j-1)*imax+i)
     end do
   end do
 endif

This is probably not the recommended way of doing things. I'm quite new
at using PETSc myself, so if anyone has a better solution, please
enlighten me! I should use PETSc data structures for everything, but I'm
trying to integrate PETSc into already existing code, so it's not that
easy to do.

Regards,
Knut Erik Teigen

> 
> As it turns out, the normal version of the system (A^t A x = A^t b)
> is still very sparse in my case. So I used MatMatMultTranspose
> and MatMultTranspose to produce the normal equations and
> solved the square system.
> 
> Using ksp_type cg and the default settings for the ilu, the solver
> converges much faster than using lsqr, which is great. However,
> MatMatMultTranspose doesn't seem to like the mpiaij matrix format.
> 
> 3) If I want to try the normal approach in multiple processes, do I
> have perform
> the product myself, or is there another function that does this for me?
> 
> Thanks in advance,
> Diego.
> 
>