[Nek5000-users] Interfacing with Arpack

[nek5000-users at lists.mcs.anl.gov]

Dear Paul,
thank you for the answer.
If my understanding is correct, since each matrix-vector product (coded 
as matrix-matrix product) is done on the local variables, I should call 
dssum after each matrix-vector operation.
For instance, if I want to apply the Helmholtz matrix to a vector u and 
store the result in v, I should do:
call axhelm(v,u,h1,h2,imsh,1)
call dssum(v,nx1,ny1,nz1)
is it correct?
If I want to apply the gradient to a vector (always for eigenvalue 
computation purposes) shoud I do something like:
       call gradm1(dudx,dudy,dudz,u)
       call dssum(dudx,nx1,ny1,nz1)
       call dssum(dudy,nx1,ny1,nz1)
       call dssum(dudz,nx1,ny1,nz1)
       bufx = glsc2(u0,dudx,n)           !bufx = u0.du/dx
or am I missing something?

Giuseppe



On 20/10/2014 12:58, nek5000-users at lists.mcs.anl.gov wrote:
> Giuseppe,
>
> The key thing to understand about Nek is that variables are _always_ stored in their local format,
> which implies you have redundant values at the interface such that they are ready for parallel operator
> evaluation.
>
> Let u_L denote such a vector with redundant values (L for Local) and let u be the corresponding
> vector without the redundantly stored data.  Nominally, there is a matrix Q that maps u to u_L:
>
> (1)       u_L = Q u
>
> and the matrix-vector product involving the mass matrix, of the form w = M u, can be expressed as:
>
> (2)       w = M u =  Q^T M_L Q u  = Q^T  M_L u_L
>
> where M_L = block_diag(M^e) , e=1,...,Nelv is the block-diagonal matrix consisting of local mass
> matrices within each element, e.   In Nek, M^e is also block diagonal, so the whole _unassembled_
> stiffness matrix, M_L, is diagonal.   It is this matrix that is stored in bm1  ("B" on Mesh 1, the velocity/temp
> mesh).
>
> Notice that, as written, the output of (2) violates our "store only local vectors" policy, so we invoke
> (1) to correct this:
>
>
> (3)    w_L = Q w = Q M u =  Q Q^T  M_L u_L.
>
>
> We refer to QQ^T as direct-stiffness summation  (dssum, in nek).   It is the operation that
> glues things together.
>
> So, suppose you wish to find a continuous function that is close to a discontinuous function
> in nek.   Let f be the continuous thing and g be the discontinuous one.   You can do this via:
>
> (4)      f_L  = Q  M^-1   Q M_L g_L.
>
> In Nek, this would be written as:
>
>            n = nx1*ny1*nz1*nelv
>
>            call col3(f,g,bm1,n)
>
>            call dssum(f,nx1,ny1,nz1)
>
>            call col2    (f,binvm1,n)
>
> Here, binvm1 is  effectively   binvm1 =  1./ (QQT*bm1)   , if we view bm1 and binvm1 as
> vectors, rather than diagonal matrices.
>
> Paul
>
>
>
>
>
> ________________________________________
> From: nek5000-users-bounces at lists.mcs.anl.gov [nek5000-users-bounces at lists.mcs.anl.gov] on behalf of nek5000-users at lists.mcs.anl.gov [nek5000-users at lists.mcs.anl.gov]
> Sent: Monday, October 20, 2014 2:03 AM
> To: nek5000-users at lists.mcs.anl.gov
> Subject: [Nek5000-users] Interfacing with Arpack
>
> Dear users and developers,
> I would like to solve a generalized eigenvalue problem with Arpack,
> based on Nek subroutines. Using reverse communication interface, I only
> need to express the action of the matrices on the vectors, but I am not
> sure on how to apply the mass matrix and the inverse mass matrix in the
> Pn-Pn formulation.
> Is it sufficient to write:
> call col2(p1,bm1,n)      !apply mass matrix to vector p1
> ...
> call invcol2(p1,bm1,n)  ! apply inverse mass matrix to vector p1
> or should I do:
> call rzero(h1,n)
> call rone(h2,n)
> call hmholtz('vel',p,p1,h1,h2,pmask,vmult,imsh,tolh,maxiter,isd)
> to apply inverse mass matrix to p1, or should I do something else?
> Also, from a related question from a couple of years ago, should I call
> col2(v,vmask,n) after each matrix-vector product?
> Thank you in advance.
> Best,
> Giuseppe
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users