Usage¶

The calculation is done to utilize functions by the following procedures.

Initialization (*_init)
Update results iteratively (*_update)
(Optional) Take the information for the restart (*_getcoef, *_getvec)
Finalization (*_finalize)

The restart calculation can be done by the following procedures.

Initialization with the information of the previous calculation (*_restart)
Update results iteratively (*_update)
(Optional) Take the information for the further restart (*_getcoef, *_getvec)
Finalization (*_finalize)

Warning

Since \(K\omega\) is not thread safe, these routine must be called from the outside of the OpenMP-parallel region. If you want to call \(K\omega\) routine in the parallel region (i.e. plan to solve different equations among threads), please use --enable-threadsafe option of configure (Experimental).

For FORTRAN, the modules can be called by

USE komega_cg_r ! Conjugate-gradient method for real vectors
USE komega_cg_c ! Conjugate-gradient method for complex vectors
USE komega_cocg ! Conjugate-orthogonal conjugate-gradient mehod
USE komega_bicg ! Biconjugate-gradient method

When we call \(K\omega\) from C/C++ codes, we should include the header file as

#include komega.h

Scalar arguments should be passed as pointers. For MPI/Hybrid parallelized routine, the above line becomes Also, the communicator argument for the routine should be transformed from the C/C++’s one to the fortran’s one as follows.

comm_f = MPI_Comm_c2f(comm_c);

Details of each routines¶

*_init¶

Set and initialize internal variables in libraries. These routines should be called first before solving the shifted equation.

Syntax

Fortran

CALL komega_cg_r_init(ndim, nl, nz, x, z, itermax, threshold, comm)
CALL komega_cg_c_init(ndim, nl, nz, x, z, itermax, threshold, comm)
CALL komega_cocg_init(ndim, nl, nz, x, z, itermax, threshold, comm)
CALL komega_bicg_init(ndim, nl, nz, x, z, itermax, threshold, comm)

C/C++

komega_cg_r_init(&ndim, &nl, &nz, x, z, &itermax, &threshold, &comm);
komega_cg_c_init(&ndim, &nl, &nz, x, z, &itermax, &threshold, &comm);
komega_cocg_init(&ndim, &nl, &nz, x, z, &itermax, &threshold, &comm);
komega_bicg_init(&ndim, &nl, &nz, x, z, &itermax, &threshold, &comm);

Parameters

INTEGER,INTENT(IN) :: ndim
The dimension of solution vectors for the linearized equation. ndim for the dimension of variables in other routine is equal to this.
INTEGER,INTENT(IN) :: nl
The dimension of projected solution vectors. nl for the dimension of variables in other routine is equal to this.
INTEGER,INTENT(IN) :: nz
The number of shifted points. nz for the dimension of variables in other routine is equal to this.
REAL(8),INTENT(OUT) :: x(nl*nz) ! (for "CG_R_init", "CG_C_init")
COMPLEX(8),INTENT(OUT) :: x(nl*nz) ! (for other cases)
The solution vector. In this procedure, 0 vector is returned.
REAL(8),INTENT(IN) :: z(nz) ! (for "CG_R_init", "CG_C_init")
COMPLEX(8),INTENT(IN) :: z(nz) ! (for other cases)
Shifted points.
INTEGER,INTENT(IN) :: itermax
The maximum iteration number for allocating arrays for the restart calculation. When itermax=0 , these arrays are not allocated, and the restart calculation described later becomes unavailable.
REAL(8),INTENT(IN) :: threshold
The threshold value for the convergence determination. When the 2-norm of the residual vector for the seed equation becomes smaller than this value, the calculation is finished.
INTEGER,INTENT(IN),OPTIONAL :: comm
Optional argument. Communicators for MPI such as MPI_COMM_WORLD. Only for MPI / Hybrid parallelization. For C compiler, just pass NULL to omit this argment.

*_restart¶

For the restart calculation, these routines are used instead of *_init. Set and initialize internal variables in libraries. These routines should be called first before solving the shifted equation.

Syntax

Fortran

CALL komega_cg_r_restart(ndim, nl, nz, x, z, itermax, threshold, status, &
&                 iter_old, v2, v12, alpha_save, beta_save, z_seed, r_l_save, comm)
CALL komega_cg_c_restart(ndim, nl, nz, x, z, itermax, threshold, status, &
&                 iter_old, v2, v12, alpha_save, beta_save, z_seed, r_l_save, comm)
CALL komega_cocg_restart(ndim, nl, nz, x, z, itermax, threshold, status, &
&                 iter_old, v2, v12, alpha_save, beta_save, z_seed, r_l_save, comm)
CALL komega_bicg_restart(ndim, nl, nz, x, z, itermax, threshold, status, &
&                 iter_old, v2, v12, v4, v14, alpha_save, beta_save, &
&                 z_seed, r_l_save, comm)

C/C++

komega_cg_r_restart(&ndim, &nl, &nz, x, z, &itermax, &threshold, status, &
&                 &iter_old, v2, v12, alpha_save, beta_save, &z_seed, r_l_save, &comm);
komega_cg_c_restart(&ndim, &nl, &nz, x, z, &itermax, &threshold, status, &
&                 &iter_old, v2, v12, alpha_save, beta_save, &z_seed, r_l_save, &comm);
komega_cocg_restart(&ndim, &nl, &nz, x, z, &itermax, &threshold, status, &
&                 &iter_old, v2, v12, alpha_save, beta_save, &z_seed, r_l_save, &comm);
komega_bicg_restart(&ndim, &nl, &nz, x, z, &itermax, &threshold, status, &
&                 &iter_old, v2, v12, v4, v14, alpha_save, beta_save, &
&                 &z_seed, r_l_save, &comm);

Parameters

INTEGER,INTENT(IN) :: ndim
INTEGER,INTENT(IN) :: nl
INTEGER,INTENT(IN) :: nz
REAL(8),INTENT(OUT) :: x(nl*nz)
REAL(8),INTENT(IN) :: z(nz) ! (for "CG_R_restart", "CG_C_restart")
COMPLEX(8),INTENT(IN) :: z(nz) ! (Other)
INTEGER,INTENT(IN) :: itermax
REAL(8),INTENT(IN) :: threshold
INTEGER,INTENT(IN),OPTIONAL :: comm
The definition is same as *_init. See the parameters in *_init.
INTEGER,INTENT(OUT) :: status(3)
The error code is returned.

First component(status(1))

If the solution is converged or a breakdown occurs, the current total number of iteration with the minus sign is returned. In other cases, this routine returns the current total number of iteration. The calculation is continuable only when status(1) is the positive value; otherwise the result is meaningless even if the calculation is continued.

Second component(status(2))

1 is returned if itermax is set as a finite value and the convergence condition is not satisfied at the itermax-th iteration. 2 is returned if \(\alpha\) diverges. 3 is returned if \(\pi_{\rm seed}\) becomes 0. In the case of COCG_restart or BiCG_restart, 4 is returned if the residual vector and the shadow residual vector are orthogonal. In other cases, 0 is returned.

Third component(status(3))

The index of the seed point is returned.
INTEGER,INTENT(IN) :: iter_old
The number of iteration for the previous calculation.
REAL(8),INTENT(IN) :: v2(ndim) ! (for "CG_R_restart")
COMPLEX(8),INTENT(IN) :: v2(ndim) ! (Other)
The residual vector at the last step for the previous calculation.
REAL(8),INTENT(IN) :: v12(ndim) ! (for "CG_R_restart")
COMPLEX(8),INTENT(IN) :: v12(ndim) ! (Other)
The residual vector at the second from the last step for the previous calculation.
REAL(8),INTENT(IN) :: alpha_save(iter_old) ! (for "CG_R_restart", "CG_C_restart")
COMPLEX(8),INTENT(IN) :: alpha_save(iter_old) ! (Other)
The parameters \(\alpha\) obtained by the previous calculation at each steps by (Bi)CG methods.
REAL(8),INTENT(IN) :: beta_save(iter_old) ! (for "CG_R_restart", "CG_C_restart")
COMPLEX(8),INTENT(IN) :: beta_save(iter_old) ! (Other)
The parameters \(\beta\) obtained by the previous calculation at each steps by (Bi)CG methods.
REAL(8),INTENT(IN) :: z_seed ! (for "CG_R_restart", "CG_C_restart")
COMPLEX(8),INTENT(IN) :: z_seed ! (Other)
The value of the seed shift for the previous calculation.
REAL(8),INTENT(IN) :: r_l_save(nl,iter_old) ! (for "CG_R_restart")
COMPLEX(8),INTENT(IN) :: r_l_save(nl,iter_old) ! (Other)
The projected residual vector at each iteration for the previous calculation.
REAL(8),INTENT(IN) :: v4(ndim) ! (for "CG_R_restart")
COMPLEX(8),INTENT(IN) :: v4(ndim) ! (Other)
The shadow residual vector at the last step for the previous calculation.
REAL(8),INTENT(IN) :: v14(ndim) ! (for "CG_R_restart")
COMPLEX(8),INTENT(IN) :: v14(ndim) ! (Other)
The shadow residual vector at the second last step for the previous calculation.

*_update¶

It is called alternately with the matrix-vector product in the loop and updates the solution.

Syntax

Fortran

CALL komega_cg_r_update(v12, v2, x, r_l, status)
CALL komega_cg_c_update(v12, v2, x, r_l, status)
CALL komega_cocg_update(v12, v2, x, r_l, status)
CALL komega_bicg_update(v12, v2, v14, v4, x, r_l, status)

C/C++

komega_cg_r_update(v12, v2, x, r_l, status);
komega_cg_c_update(v12, v2, x, r_l, status);
komega_cocg_update(v12, v2, x, r_l, status);
komega_bicg_update(v12, v2, v14, v4, x, r_l, status);

Parameters

REAL(8),INTENT(INOUT) :: v12(ndim) ! (for "CG_R_update")
COMPLEX(8),INTENT(INOUT) :: v12(ndim) ! (Other)
The product of the residual vector (v2) and the matrix. This routine returns the 2-norm of the updated residual vector as a first element of this array. This returned value is used, for examples, for printing the convergence profile.
REAL(8),INTENT(INOUT) :: v2(ndim) ! (for "CG_R_update")
COMPLEX(8),INTENT(INOUT) :: v2(ndim) ! (Other)
The residual vector is input and the updated residual vector is output.
REAL(8),INTENT(IN) :: v14(ndim) ! (for "CG_R_update")
COMPLEX(8),INTENT(IN) :: v14(ndim) ! (Other)
The product of the shadow residual vector (v4) and the matrix is input.
REAL(8),INTENT(INOUT) :: v4(ndim) ! (for "CG_R_update")
COMPLEX(8),INTENT(INOUT) :: v4(ndim) ! (Other)
The shadow residual vector is input and the updated vector is output.
INTEGER,INTENT(OUT) :: status(3)
The error code is returned.

First component (status(1))

If the solution is converged or a breakdown occurs, the current total number of iteration with the minus sign is returned. In other cases, this routine returns the current total number of iteration. The calculation is continuable only when status(1) is the positive value; otherwise the result is meaningless even if the calculation is continued.

Second component (status(2))

1 is returned if itermax is set as a finite value in the *_init routine and the convergence condition is not satisfied at the itermax-th iteration. 2 is returned if \(\alpha\) diverges. 3 is returned if \(\pi_{\rm seed}\) becomes 0. In the case of COCG_update or BiCG_update, 4 is returned if the residual vector and the shadow residual vector are orthogonal. In other cases, 0 is returned.

Third component (status(3))

The index of the seed point is returned.

*_getcoef¶

Get the coefficients used in the restart calculation. To call these routines, itermax in *_init routine must not be 0 .

The total number of iteration (iter_old) used in this routine is computed by using status which is an output of *_update as follows:

iter_old = ABS(status(1))

Syntax

Fortran

CALL komega_cg_r_getcoef(alpha_save, beta_save, z_seed, r_l_save)
CALL komega_cg_c_getcoef(alpha_save, beta_save, z_seed, r_l_save)
CALL komega_cocg_getcoef(alpha_save, beta_save, z_seed, r_l_save)
CALL komega_bicg_getcoef(alpha_save, beta_save, z_seed, r_l_save)

C/C++

komega_cg_r_getcoef(alpha_save, beta_save, &z_seed, r_l_save);
komega_cg_c_getcoef(alpha_save, beta_save, &z_seed, r_l_save);
komega_cocg_getcoef(alpha_save, beta_save, &z_seed, r_l_save);
komega_bicg_getcoef(alpha_save, beta_save, &z_seed, r_l_save);

Parameters

REAL(8),INTENT(OUT) :: alpha_save(iter_old) ! (for "CG_R_restart", "CG_C_restart")
COMPLEX(8),INTENT(OUT) :: alpha_save(iter_old) ! (Other)

The parameters \(\alpha\) of the (Bi)CG method at each iteration.

REAL(8),INTENT(OUT) :: beta_save(iter_old) ! (for "CG_R_restart", "CG_C_restart")
COMPLEX(8),INTENT(OUT) :: beta_save(iter_old) ! (Other)

The parameters \(\beta\) of the (Bi)CG method at each iteration.

REAL(8),INTENT(OUT) :: z_seed ! (for "CG_R_restart", "CG_C_restart")
COMPLEX(8),INTENT(OUT) :: z_seed ! (Other)

Seed shift.

REAL(8),INTENT(IN) :: r_l_save(nl,iter_old) ! ("CG_R_restart")
COMPLEX(8),INTENT(IN) :: r_l_save(nl,iter_old) ! (Other)

The projected residual vectors at each iteration.

*_getvec¶

Get the residual vectors to use the restart calculation. To call these routines, itermax in the *_init routine must not be 0.

Syntax

Fortran

CALL komega_cg_r_getvec(r_old)
CALL komega_cg_c_getvec(r_old)
CALL komega_cocg_getvec(r_old)
CALL komega_bicg_getvec(r_old, r_tilde_old)

C/C++

komega_cg_r_getvec(r_old);
komega_cg_c_getvec(r_old);
komega_cocg_getvec(r_old);
komega_bicg_getvec(r_old, r_tilde_old);

Parameters

REAL(8),INTENT(OUT) :: r_old(ndim) ! (for "CG_R_getvec")
COMPLEX(8),INTENT(OUT) :: r_old(ndim) ! (Other)
The residual vector at the second last step in the previous calculation.
COMPLEX(8),INTENT(OUT) :: r_tilde_old(ndim)
The shadow residual vector at the second last step in the previous calculation.

*_getresidual¶

Get the values of 2-norm of the residual vector at each shift points. These routines can be called from anywhere between *_init and *_finalize . These routines do not affect the calculation results.

Syntax

Fortran

CALL komega_cg_r_getresidual(res)
CALL komega_cg_c_getresidual(res)
CALL komega_cocg_getresidual(res)
CALL komega_bicg_getresidual(res)

C/C++

komega_cg_r_getresidual(res);
komega_cg_c_getresidual(res);
komega_cocg_getresidual(res);
komega_bicg_getresidual(res);

Parameters

COMPLEX(8),INTENT(OUT) :: res(nz)
The values of 2-norm of the residual vector at each shift points are returned.

*_finalize¶

Release memories of the arrays stored in the library.

Syntax

Fortran

CALL komega_cg_r_finalize()
CALL komega_cg_c_finalize()
CALL komega_cocg_finalize()
CALL komega_bicg_finalize()

C/C++

komega_cg_r_finalize();
komega_cg_c_finalize();
komega_cocg_finalize();
komega_bicg_finalize();

Sample codes for using shifted BiCG library¶

As a typical example, the usage of shifted BiCG library is shown below.

PROGRAM my_prog
  !
  USE komega_bicg, ONLY : komega_bicg_init, komega_bicg_restart, &
  &                       komega_bicg_update, komega_bicg_getcoef, &
  &                       komega_bicg_getvec, komega_bicg_finalize
  USE solve_cc_routines, ONLY : input_size, input_restart, &
  &                             projection, &
  &                             hamiltonian_prod, generate_system, &
  &                             output_restart, output_result
  !
  IMPLICIT NONE
  !
  INTEGER,SAVE :: &
  & ndim,    & ! Size of Hilvert space
  & nz,      & ! Number of frequencies
  & nl,      & ! Number of Left vector
  & itermax, & ! Max. number of iteraction
  & iter_old   ! Number of iteraction of previous run
  !
  REAL(8),SAVE :: &
  & threshold ! Convergence Threshold
  !
  COMPLEX(8),SAVE :: &
  & z_seed ! Seed frequency
  !
  COMPLEX(8),ALLOCATABLE,SAVE :: &
  & z(:)         ! (nz): Frequency
  !
  COMPLEX(8),ALLOCATABLE,SAVE :: &
  & ham(:,:), &
  & rhs(:), &
  & v12(:), v2(:), & ! (ndim): Working vector
  & v14(:), v4(:), & ! (ndim): Working vector
  & r_l(:), & ! (nl) : Projeccted residual vector
  & x(:,:) ! (nl,nz) : Projected result
  !
  ! Variables for Restart
  !
  COMPLEX(8),ALLOCATABLE,SAVE :: &
  & alpha(:), beta(:) ! (iter_old)
  !
  COMPLEX(8),ALLOCATABLE,SAVE :: &
  & r_l_save(:,:) ! (nl,iter_old) Projected residual vectors
  !
  ! Variables for Restart
  !
  INTEGER :: &
  & iter,    & ! Counter for Iteration
  & status(3)
  !
  LOGICAL :: &
  & restart_in, & ! If .TRUE., sestart from the previous result
  & restart_out   ! If .TRUE., save datas for the next run
  !
  ! Input Size of vectors, numerical conditions
  !
  CALL input_size(ndim,nl,nz)
  CALL input_condition(itermax,threshold,restart_in,restart_out)
  !
  ALLOCATE(v12(ndim), v2(ndim), v14(ndim), v4(ndim), r_l(nl), &
  &        x(nl,nz), z(nz), ham(ndim,ndim), rhs(ndim))
  !
  CALL generate_system(ndim, ham, rhs, z)
  !
  WRITE(*,*)
  WRITE(*,*) "#####  CG Initialization  #####"
  WRITE(*,*)
  !
  IF(restart_in) THEN
    !
    CALL input_restart(iter_old, zseed, alpha, beta, r_l_save)
    !
    IF(restart_out) THEN
       CALL komega_bicg_restart( &
       &    ndim, nl, nz, x, z, itermax, threshold, &
       &    status, iter_old, v2, v12, v4, v14, alpha, &
       &    beta, z_seed, r_l_save)
    ELSE
       CALL komega_bicg_restart( &
       &    ndim, nl, nz, x, z, 0, threshold, &
       &    status, iter_old, v2, v12, v4, v14, alpha, &
       &    beta, z_seed, r_l_save)
    END IF
    !
    ! These vectors were saved in BiCG routine
    !
    DEALLOCATE(alpha, beta, r_l_save)
    !
    IF(status(1) /= 0) GOTO 10
    !
  ELSE
     !
     ! Generate Right Hand Side Vector
     !
     v2(1:ndim) = rhs(1:ndim)
     v4(1:ndim) = CONJG(v2(1:ndim))
     !v4(1:ndim) = v2(1:ndim)
     !
     IF(restart_out) THEN
        CALL komega_bicg_init(ndim, nl, nz, x, z, termax, threshold)
     ELSE
        CALL komega_bicg_init(ndim, nl, nz, x, z, 0, threshold)
     END IF
     !
  END IF
  !
  ! BiCG Loop
  !
  WRITE(*,*)
  WRITE(*,*) "#####  CG Iteration  #####"
  WRITE(*,*)
  !
  DO iter = 1, itermax
     !
     ! Projection of Residual vector into the space
     ! spaned by left vectors
     !
     r_l(1:nl) = projection(v2(1:nl))
     !
     ! Matrix-vector product
     !
     CALL hamiltonian_prod(Ham, v2, v12)
     CALL hamiltonian_prod(Ham, v4, v14)
     !
     ! Update result x with BiCG
     !
     CALL komega_bicg_update(v12, v2, v14, v4, x, r_l, status)
     !
     WRITE(*,'(a,i,a,3i,a,e15.5)') "lopp : ", iter, &
     &                             ", status : ", status(1:3), &
     &                             ", Res. : ", DBLE(v12(1))
     IF(status(1) < 0) EXIT
     !
  END DO
  !
  IF(status(2) == 0) THEN
     WRITE(*,*) "  Converged in iteration ", ABS(status(1))
  ELSE IF(status(2) == 1) THEN
     WRITE(*,*) "  Not Converged in iteration ", ABS(status(1))
  ELSE IF(status(2) == 2) THEN
     WRITE(*,*) "  Alpha becomes infinity", ABS(status(1))
  ELSE IF(status(2) == 3) THEN
     WRITE(*,*) "  Pi_seed becomes zero", ABS(status(1))
  ELSE IF(status(2) == 4) THEN
  WRITE(*,*) "  Residual & Shadow residual are orthogonal", &
  &          ABS(status(1))
  END IF
  !
  ! Total number of iteration
  !
  iter_old = ABS(status(1))
  !
  ! Get these vectors for restart in the Next run
  !
  IF(restart_out) THEN
     !
     ALLOCATE(alpha(iter_old), beta(iter_old), r_l_save(nl,iter_old))
     !
     CALL komega_bicg_getcoef(alpha, beta, z_seed, r_l_save)
     CALL komega_bicg_getvec(v12,v14)
     !
     CALL output_restart(iter_old, z_seed, alpha, beta, &
     &                   r_l_save, v12, v14)
     !
     DEALLOCATE(alpha, beta, r_l_save)
     !
  END IF
  !
10 CONTINUE
  !
  ! Deallocate all intrinsic vectors
  !
  CALL komega_bicg_finalize()
  !
  ! Output to a file
  !
  CALL output_result(nl, nz, z, x, r_l)
  !
  DEALLOCATE(v12, v2, v14, v4, r_l, x, z)
  !
  WRITE(*,*)
  WRITE(*,*) "#####  Done  #####"
  WRITE(*,*)
  !
END PROGRAM my_prog

Usage¶

Details of each routines¶

*_init¶

*_restart¶

*_update¶

*_getcoef¶

*_getvec¶

*_getresidual¶

*_finalize¶

Sample codes for using shifted BiCG library¶

Table Of Contents

Previous topic

Next topic

This Page

Navigation

Usage¶

Details of each routines¶

*_init¶

*_restart¶

*_update¶

*_getcoef¶

*_getvec¶

*_getresidual¶

*_finalize¶

Sample codes for using shifted BiCG library¶

Table Of Contents

Previous topic

Next topic

This Page

Quick search

Navigation