matrix free implementation issue

[felipeZ]

In the current implementation of the Matrix-free version:

• The projected matrix is updated by block but the QR factorization could change the sign of the vectors making a block-update unfeasible.
• It was assumed that matrix multiplication of a diagonal matrix is commutative that is wrong for nonsymmetric matrices.
• The correction matrix is computed with the wrong diagonal matrix because for each eigenvalue/eigenvector pair a diagonal matrix must be formed containing a single eigenvalue in the diagonal.
• The calls to the function that computes the matrix ( the matrix-free component) are duplicated in several locations

Todo:

• [x] Remove the block update
• [x] Fix multiplication matrix
• [x] Fix diagonal matrix
• [x] Optimize the number of calls to the function that computes the target matrix.

[felipeZ]

@lopeztarifa I have done a couple of optimizations and fix the issues mention above. Do you think that there are further improvements to do?

[lopeztarifa]

Great work!

Now the compute_DPR_free is much more elegant and efficient.

Just a final remark, we discussed the addition of a third argument to the generate_preconditioner function so you can change the first dimension of precond. It is not there, right? Might be it was to specif to the Champ case ... it is Friday and I am off guard. This is what I mean:

Free: _V = generate_preconditioner( d(1:dim_sub), dimsub, nparm) ! Initial orthonormal basis Desnse: _V = generate_preconditioner(d(1:dim_sub), dimsub, size(mtx,1))

where:

  _function generate_preconditioner(diag, dim_sub, dim_base) result(precond)
    !> \brief generates a diagonal preconditioner for `matrix`.
    !> \return diagonal matrix

    ! input variable
    real(dp), dimension(:), intent(inout) :: diag
    integer, intent(in) :: dim_sub, dim_base

    ! local variables
    real(dp), dimension(dim_base, dim_sub) :: precond                                                                                                                                                              
    integer, dimension(size(diag)) :: keys
    integer :: i, k

    ! sort diagonal
    keys = lapack_sort('I', diag)
    ! Fill matrix with zeros
    precond = 0.0_dp

    ! Add one depending on the order of the matrix diagonal
    do i=1, dim_sub
       k = search_key(keys, i)
       precond(k, i) = 1.d0
    end do

  end function generate_preconditioner_

[felipeZ]

The generate_preconditioner subroutine is already called in this line for the free implementation and define here.