module LinearSolverClass use, intrinsic :: iso_fortran_env use omp_lib use IOClass use TimeClass use MathClass use ArrayClass use SparseClass use RandomClass !use MPIClass implicit none interface BiCGSTAB module procedure bicgstab_real32, bicgstab_real64, bicgstab_complex64 end interface interface GPBiCG module procedure GPBiCG_real32, GPBiCG_real64, GPBiCG_complex64 end interface interface Gauss_Jordan_PV module procedure Gauss_Jordan_PV_real32, Gauss_Jordan_PV_real64, Gauss_Jordan_PV_complex64 end interface interface StochasticGradientDescent module procedure StochasticGradientDescent_scalar_function end interface type :: LinearSolver_ ! non-Element-by-element real(real64), allocatable :: a(:, :) real(real64), allocatable :: b(:) real(real64), allocatable :: x(:) ! Element-by-element real(real64), allocatable :: a_e(:, :, :) real(real64), allocatable :: b_e(:, :) real(real64), allocatable :: x_e(:, :) ! For Sparse ! COO format real(real64), allocatable :: val(:) integer(int32), allocatable :: index_I(:) integer(int32), allocatable :: index_J(:) integer(int32), allocatable :: row_domain_id(:) integer(int32), allocatable :: column_domain_id(:) ! CRS format real(real64), allocatable :: CRS_val(:) integer(int64), allocatable :: CRS_index_I(:) integer(int32), allocatable :: CRS_index_J(:) integer(int32), allocatable :: CRS_row_domain_id(:) integer(int32), allocatable :: CRS_column_domain_id(:) integer(int32), allocatable :: b_Index_J(:) integer(int32), allocatable :: b_Domain_ID(:) logical, allocatable :: Locked(:) ! info integer(int32), allocatable :: NumberOfNode(:) integer(int32) :: DOF = 1 logical :: debug = .false. ! integer(int32), allocatable :: connectivity(:, :) integer(int32) :: itrmax = 1000000 integer(int32) :: currentID = 1 integer(int32) :: b_currentID = 1 real(real64) :: er0 = dble(1.0e-10) logical :: ReadyForFix = .false. contains procedure, public :: init => initLinearSolver procedure, public :: set => setLinearSolver procedure, public :: assemble => assembleLinearSolver procedure, pass :: importLinearSolver procedure, pass :: importCRS_LinearSolver generic, public :: import => importLinearSolver, importCRS_LinearSolver procedure, public :: fix => fixLinearSolver procedure, public :: solve => solveLinearSolver procedure, public :: show => showLinearSolver procedure, public :: globalMatrix => globalMatrixLinearSolver procedure, public :: globalVector => globalVectorLinearSolver procedure, public :: convertCOOtoCRS => convertCOOtoCRSLinearSolver procedure, public :: matmulCRS => matmulCRSLinearSolver procedure, public :: matmulCOO => matmulCOOLinearSolver procedure, public :: prepareFix => prepareFixLinearSolver procedure, public :: getCOOFormat => getCOOFormatLinearSolver procedure, public :: exportAsCOO => exportAsCOOLinearSolver procedure, public :: exportRHS => exportRHSLinearSolver end type interface crs_matvec module procedure :: crs_matvec_generic, crs_matvec_for_CRStype end interface crs_matvec interface sub_crs_matvec module procedure :: sub_crs_matvec_generic, sub_crs_matvec_for_CRStype end interface sub_crs_matvec interface crs_matmul module procedure :: crs_matmul_generic, crs_matmul_for_CRStype end interface crs_matmul contains !==================================================================================== subroutine initLinearSolver(obj, NumberOfNode, DOF) class(LinearSolver_), intent(inout) :: obj integer(int32), optional, intent(in) :: NumberOfNode(:), DOF integer(int32) :: i, j, k, n, num_total_unk, node_count obj%ReadyForFix = .false. ! non-Element-by-element if (allocated(obj%a)) deallocate (obj%a) if (allocated(obj%b)) deallocate (obj%b) if (allocated(obj%x)) deallocate (obj%x) ! Element-by-element if (allocated(obj%a_e)) deallocate (obj%a_e) if (allocated(obj%b_e)) deallocate (obj%b_e) if (allocated(obj%x_e)) deallocate (obj%x_e) if (allocated(obj%val)) deallocate (obj%val) if (allocated(obj%index_I)) deallocate (obj%index_I) if (allocated(obj%index_J)) deallocate (obj%index_J) if (allocated(obj%row_Domain_ID)) deallocate (obj%row_Domain_ID) if (allocated(obj%column_Domain_ID)) deallocate (obj%column_Domain_ID) if (allocated(obj%b_Index_J)) deallocate (obj%b_Index_J) if (allocated(obj%b_Domain_ID)) deallocate (obj%b_Domain_ID) if (allocated(obj%Locked)) deallocate (obj%Locked) if (allocated(obj%connectivity)) deallocate (obj%connectivity) n = input(default=1, option=DOF) ! Number of node of n th domains is NumberOfNode(n) if (present(NumberOfNode)) then num_total_unk = sum(NumberOfNode)*n allocate (obj%b_Index_J(num_total_unk)) allocate (obj%b_Domain_ID(num_total_unk)) allocate (obj%Locked(num_total_unk)) allocate (obj%b(num_total_unk)) obj%NumberOfNode = NumberOfNode obj%DOF = DOF obj%b(:) = 0.0d0 num_total_unk = 0 do i = 1, size(NumberOfNode) node_count = 0 do j = 1, NumberOfNode(i) do k = 1, n num_total_unk = num_total_unk + 1 node_count = node_count + 1 obj%b_Domain_ID(num_total_unk) = i obj%Locked(num_total_unk) = .false. obj%b_Index_J(num_total_unk) = node_count end do end do end do end if obj%itrmax = 1000000 obj%currentID = 1 obj%b_currentID = 1 obj%er0 = dble(1.0e-08) end subroutine !==================================================================================== !==================================================================================== recursive subroutine assembleLinearSolver(obj, connectivity, DOF, eMatrix, eVector, DomainIDs) class(LinearSolver_), intent(inout) :: obj integer(int32), intent(in) :: connectivity(:) ! connectivity matrix !(global_node_id#1, global_node_id#2, global_node_id#3, . ) integer(int32), intent(in) :: DOF ! degree of freedom integer(int32), optional, intent(in) :: DomainIDs(:) ! DomainID real(real64), optional, intent(in) :: eMatrix(:, :) ! elemental matrix real(real64), optional, intent(in) :: eVector(:) ! elemental Vector integer(int32) :: i, j, k, l, m, node_id1, node_id2, domain_ID1, domain_ID2 integer(int32), allocatable :: domID(:) if (present(eMatrix)) then if (present(DomainIDs)) then do j = 1, size(connectivity, 1) do k = 1, size(connectivity, 1) do l = 1, DOF do m = 1, DOF node_id1 = connectivity(j) node_id2 = connectivity(k) if (j < 1 .or. k < 1) then print *, "ERROR :: Assemble solver j<1 .or. k<1" stop end if if (j > size(DomainIDs) .or. k > size(DomainIDs)) then print *, j, size(DomainIDs), k print *, DOmainIDs print *, "ERROR :: Assemble solver j >= size(DomainIDs) .or.k >= size(DomainIDs)" stop end if domain_ID1 = DomainIDs(j) domain_ID2 = DomainIDs(k) call obj%set( & low=DOF*(node_id1 - 1) + l, & column=DOF*(node_id2 - 1) + m, & entryvalue=eMatrix(DOF*(j - 1) + l, DOF*(k - 1) + m), & row_DomainID=Domain_ID1, & column_DomainID=Domain_ID2) end do end do end do end do else do j = 1, size(DomainIDs, 1) do k = 1, size(DomainIDs, 1) do l = 1, DOF do m = 1, DOF node_id1 = connectivity(j) node_id2 = connectivity(k) call obj%set( & low=DOF*(node_id1 - 1) + l, & column=DOF*(node_id2 - 1) + m, & entryvalue=eMatrix(DOF*(j - 1) + l, DOF*(k - 1) + m), & row_DomainID=1, & column_DomainID=1) end do end do end do end do end if end if if (present(eVector)) then if (present(DomainIDs)) then do j = 1, size(connectivity) do l = 1, DOF node_id1 = connectivity(j) domain_ID1 = DomainIDs(j) call obj%set( & low=DOF*(node_id1 - 1) + l, & entryvalue=eVector(DOF*(j - 1) + l), & row_DomainID=Domain_ID1) end do end do else do j = 1, size(connectivity) do l = 1, DOF node_id1 = connectivity(j) call obj%set( & low=DOF*(node_id1 - 1) + l, & entryvalue=eVector(DOF*(j - 1) + l), & row_DomainID=1) end do end do end if end if end subroutine !==================================================================================== !==================================================================================== recursive subroutine fixLinearSolver(obj, nodeid, entryvalue, entryID, DOF, row_DomainID, & nodeids, entryvalues, debug) class(LinearSolver_), intent(inout) :: obj logical, optional, intent(in) :: debug integer(int32), intent(in) :: nodeid integer(int32), optional, intent(in) :: entryID, DOF, row_DomainID, nodeids(:), entryvalues(:) real(real64), intent(in) :: entryvalue integer(int32), allocatable :: Index_I(:), Index_J(:), NumNodeBeforeDomainID(:) integer(int32) :: i, j, n, offset, m type(Time_) :: time ! [CASE1] :: single-domain & dense matrix ! [CASE2] :: single-domain & Sparse matrix ! [CASE3] :: multi-domain & dense matrix ! [CASE4] :: multi-domain & Sparse matrix ! [CASE0] :: none of above ! [CASE1][CASE0] if (.not. allocated(obj%val)) then if (allocated(obj%a) .and. allocated(obj%b)) then ! [CASE1] ! it only has obj%a and obj%b !x(nodeid) = entryvalue n = size(obj%b) obj%b(:) = obj%b(:) - obj%a(:, nodeid)*entryvalue obj%a(:, nodeid) = 0.0d0 obj%a(nodeid, :) = 0.0d0 obj%a(nodeid, nodeid) = 1.0d0 obj%b(nodeid) = entryvalue return else ![CASE0] print *, "[ERROR] LinearSolver >> [CASE0] :: No Ax=b is set." end if end if ! Below, cases [CASE2] or [CASE4] if (.not. present(row_DomainID)) then ![CASE2] >> set : row_DomainID = 1 >> [Case4] call obj%fix(nodeid=nodeid, entryvalue=entryvalue, entryID=entryID, DOF=DOF, & row_DomainID=1, debug=debug) end if if (obj%debug) then print *, "[fixLinearSolver] >> ReadyForFix" end if if (.not. obj%ReadyForFix) then call obj%prepareFix() end if if (obj%debug) then print *, "[Done] >> ReadyForFix" end if if (present(debug)) then if (debug) then call time%start() print *, "fixLinearSolver >> [0] started!" end if end if if (present(DOF)) then if (.not. present(entryID)) then print *, "ERROR :: fixLinearSolver >> argument [DOF] should be called with [entryID]" print *, "e.g. x-entry of nodeid=10 in terms of 3d(x,y,z) space is >> " print *, "nodeid =10, entryID=1, DOF=3" stop end if n = (nodeid - 1)*DOF + entryID if (present(row_DomainID)) then call obj%fix(nodeid=n, entryvalue=entryvalue, row_DomainID=row_DomainID, debug=debug) else call obj%fix(nodeid=n, entryvalue=entryvalue, row_DomainID=1, debug=debug) end if return end if if (.not. allocated(obj%row_domain_id) .and. .not. present(row_DomainID)) then ! only for COO-format if (.not. allocated(obj%val) .or. .not. allocated(obj%b)) then print *, "ERROR >> fixLinearSolver .not. allocated(val) " stop end if if (present(debug)) then if (debug) then print *, "fixLinearSolver >> [1] fix started!" end if end if do i = 1, size(obj%val) if (obj%index_J(i) == nodeid) then if (.not. obj%Locked(obj%index_I(i))) then obj%b(obj%index_I(i)) = obj%b(obj%index_I(i)) - obj%val(i)*entryvalue end if obj%val(i) = 0.0d0 end if end do do i = 1, size(obj%index_I) if (obj%index_I(i) == nodeid) then if (obj%index_J(i) == nodeid) then obj%val(i) = 1.0d0 else obj%val(i) = 0.0d0 end if if (.not. obj%Locked(obj%index_I(i))) then obj%b(obj%index_I(i)) = entryvalue obj%Locked(obj%index_I(i)) = .true. end if !print *, "Locked ",obj%index_I(i),"by",entryvalue end if if (obj%index_I(i) == nodeid) then if (obj%index_J(i) == nodeid) then obj%val(i) = 1.0d0 end if end if end do if (present(debug)) then if (debug) then print *, "fixLinearSolver >> [ok] done!" end if end if return elseif (allocated(obj%row_domain_id) .and. present(row_DomainID)) then ! only for COO-format if (present(debug)) then if (debug) then print *, "fixLinearSolver >> [1] Multi-domain started!" call time%show() end if end if ! fix b vector n = row_DomainID if (n == 1) then offset = 0 else offset = sum(obj%NumberOfNode(1:n - 1))*obj%DOF end if if (.not. obj%Locked(offset + nodeid)) then obj%b(offset + nodeid) = entryvalue obj%Locked(offset + nodeid) = .true. !print *, "Locked ",(offset + nodeid),"by",entryvalue end if if (present(debug)) then if (debug) then print *, "fixLinearSolver >> [1-2] count offset number" call time%show() end if end if allocate (NumNodeBeforeDomainID(size(obj%row_Domain_ID, 1))) NumNodeBeforeDomainID(1) = 0 do m = 2, size(obj%NumberOfNode) NumNodeBeforeDomainID(m) = sum(obj%NumberOfNode(1:m - 1)) end do if (.not. allocated(obj%val) .or. .not. allocated(obj%b)) then print *, "ERROR >> fixLinearSolver .not. allocated(val) " stop end if !print *, "obj%b",obj%b if (present(debug)) then if (debug) then print *, "fixLinearSolver >> [2] Updating b-vector" call time%show() end if end if ! update b-vector (Right-hand side vector) do i = 1, size(obj%val) if (obj%val(i) == 0.0d0) then cycle end if if (obj%column_Domain_ID(i) /= row_DomainID) then cycle end if if (obj%index_J(i) == nodeid) then n = obj%row_Domain_ID(i) offset = NumNodeBeforeDomainID(n)*obj%DOF n = obj%Index_I(i) if (.not. obj%Locked(offset + n)) then if (size(obj%NumberOfNode) == 1) then obj%b(n) = obj%b(n) - obj%val(i)*entryvalue else obj%b(offset + n) = obj%b(offset + n) - obj%val(i)*entryvalue end if end if obj%val(i) = 0.0d0 else cycle end if end do if (present(debug)) then if (debug) then print *, "fixLinearSolver >> [3] Updated b-vector" call time%show() end if end if do i = 1, size(obj%index_I) ! for all queries of A matrix if (obj%index_I(i) == nodeid .and. obj%row_Domain_ID(i) == row_DomainID) then obj%val(i) = 0.0d0 if (obj%index_J(i) == nodeid .and. obj%column_Domain_ID(i) == row_DomainID) then obj%val(i) = 1.0d0 end if end if end do if (present(debug)) then if (debug) then print *, "fixLinearSolver >> [ok] Done!" call time%show() end if end if else print *, "ERROR :: fixLinearSolver >> allocated(obj%row_domain_id) /= present(row_DomainID)" print *, allocated(obj%row_domain_id), present(row_DomainID) stop end if end subroutine !==================================================================================== !==================================================================================== recursive subroutine setLinearSolver(obj, low, column, entryvalue, init, row_DomainID, column_DomainID) class(LinearSolver_), intent(inout) :: obj integer(int32), optional, intent(in) :: low, column, row_DomainID, column_DomainID real(real64), optional, intent(in) :: entryvalue logical, optional, intent(in) :: init integer(int32) :: i, row_DomID, column_DomID, row_offset, column_offset, j, k, find_num, max_thread row_DomID = input(default=1, option=row_DomainID) column_DomID = input(default=1, option=column_DomainID) if (allocated(obj%NumberOfNode)) then if (present(row_DomainID)) then if (row_DomainID == 1) then row_offset = 0 else row_offset = sum(obj%NumberOfNode(1:row_DomainID - 1)) end if end if if (present(column_DomainID)) then if (column_DomainID == 1) then column_offset = 0 else column_offset = sum(obj%NumberOfNode(1:column_DomainID - 1)) end if end if end if if (present(init)) then if (init .eqv. .true.) then if (allocated(obj%val)) then obj%val(:) = 0.0d0 end if if (allocated(obj%b)) then obj%b(:) = 0.0d0 end if obj%currentID = 1 end if end if if (present(low) .and. present(column)) then if (.not. allocated(obj%val)) then allocate (obj%val(1)) obj%val(1) = input(default=0.0d0, option=entryvalue) allocate (obj%index_I(1)) obj%index_I(1) = low allocate (obj%index_J(1)) obj%index_J(1) = column allocate (obj%row_Domain_ID(1)) allocate (obj%column_Domain_ID(1)) obj%row_Domain_ID(obj%currentID) = row_DomID obj%column_Domain_ID(obj%currentID) = column_DomID ! DomainIDは，rowとcolumnの両方必要！！！ ! DomainIDは，rowとcolumnの両方必要！！！ ! DomainIDは，rowとcolumnの両方必要！！！ ! DomainIDは，rowとcolumnの両方必要！！！ ! DomainIDは，rowとcolumnの両方必要！！！ ! DomainIDは，rowとcolumnの両方必要！！！ ! DomainIDは，rowとcolumnの両方必要！！！ return end if ! if already exists, add. ! find_num = 0 ! !$omp parallel ! !$omp do reduction(+:find_num) ! do i=1,size(obj%index_I) ! if(obj%index_i(i)==0 ) cycle ! if(obj%row_domain_id(i) == row_DomainID .and. obj%column_domain_id(i) == column_DomainID )then ! if(obj%index_I(i) == low )then ! if(obj%index_J(i) == column )then ! obj%val(i) = obj%val(i) + entryvalue ! find_num = 1 ! endif ! endif ! endif ! enddo ! !$omp end do ! !$omp end parallel ! if(find_num>=1)then ! return ! endif if (obj%currentID + 1 > size(obj%val)) then call extendArray(obj%val, 0.0d0, size(obj%val)) call extendArray(obj%index_I, 0, size(obj%index_I)) call extendArray(obj%index_J, 0, size(obj%index_J)) call extendArray(obj%row_Domain_ID, 0, size(obj%row_Domain_ID)) call extendArray(obj%column_Domain_ID, 0, size(obj%column_Domain_ID)) end if obj%currentID = obj%currentID + 1 obj%val(obj%currentID) = entryvalue obj%index_I(obj%currentID) = low obj%index_J(obj%currentID) = column if (present(row_DomainID)) then obj%row_Domain_ID(obj%currentID) = row_DomainID else obj%row_Domain_ID(obj%currentID) = 1 end if if (present(column_DomainID)) then obj%column_Domain_ID(obj%currentID) = column_DomainID else obj%column_Domain_ID(obj%currentID) = 1 end if return elseif (present(low) .and. .not. present(column)) then ! for right-hand side vector if (present(row_DomainID)) then ! multi-domain problem if (.not. allocated(obj%b)) then print *, "ERROR :: setLinearSolver >> for multi-domain problem, please call %init method" print *, "with % init( NumberOfNode , DOF )" stop else ! Right-hand side vector ! Extend one-by-one if (row_DomainID == 1) then row_offset = 0 else row_offset = sum(obj%NumberOfNode(1:row_DomainID - 1))*obj%DOF end if obj%b(row_offset + low) = obj%b(row_offset + low) + entryvalue end if return else ! single-domain problem if (.not. allocated(obj%b)) then allocate (obj%b(low)) obj%b(low) = input(default=0.0d0, option=entryvalue) obj%CurrentID = low return else ! Right-hand side vector if (low > size(obj%b)) then if (obj%currentID < size(obj%val)) then obj%b(low) = entryvalue else call extendArray(obj%b, 0.0d0, low - size(obj%b)) obj%b(low) = obj%b(low) + entryvalue end if end if end if end if else return end if end subroutine !==================================================================================== !==================================================================================== subroutine importLinearSolver(obj, a, x, b, a_e, b_e, x_e, connectivity, val, index_I, index_J) class(LinearSolver_), intent(inout) :: obj real(8), optional, intent(in) :: a(:, :), b(:), x(:), a_e(:, :, :), b_e(:, :), x_e(:, :) real(8), optional, intent(in) :: val(:) integer(int32), optional, intent(in) :: index_I(:), index_J(:) integer(int32), optional, intent(in) :: connectivity(:, :) integer(int32) :: k, l, m if (present(val)) then if (.not. allocated(obj%val)) then allocate (obj%val(size(val))) obj%val = val end if end if if (present(index_i)) then if (.not. allocated(obj%index_i)) then allocate (obj%index_i(size(index_i))) obj%index_i = index_i end if end if if (present(index_j)) then if (.not. allocated(obj%index_j)) then allocate (obj%index_j(size(index_j))) obj%index_j = index_j end if end if ! in case of non element-by-element if (present(a)) then ! Set Ax=b k = size(a, 1) l = size(a, 2) if (.not. allocated(obj%a)) then allocate (obj%a(k, l)) elseif (size(obj%a, 1) /= k .or. size(obj%a, 2) /= l) then deallocate (obj%a) allocate (obj%a(k, l)) end if obj%a(:, :) = a(:, :) end if if (present(b)) then k = size(b, 1) if (.not. allocated(obj%b)) then allocate (obj%b(k)) elseif (size(obj%b, 1) /= k) then deallocate (obj%b) allocate (obj%b(k)) end if obj%b(:) = b(:) end if if (present(x)) then k = size(x, 1) if (.not. allocated(obj%x)) then allocate (obj%x(k)) elseif (size(obj%x, 1) /= k) then deallocate (obj%x) allocate (obj%x(k)) end if obj%x(:) = x(:) end if if (present(a_e)) then ! Set A_e x_e =b_e k = size(a_e, 1) l = size(a_e, 2) m = size(a_e, 3) if (.not. allocated(obj%a_e)) then allocate (obj%a_e(k, l, m)) end if obj%a_e(:, :, :) = a_e(:, :, :) end if if (present(b_e)) then ! Set Ax=b k = size(b_e, 1) l = size(b_e, 2) if (.not. allocated(obj%b_e)) then allocate (obj%b_e(k, l)) elseif (size(obj%b_e, 1) /= k .or. size(obj%b_e, 2) /= l) then deallocate (obj%b_e) allocate (obj%b_e(k, l)) end if obj%b_e(:, :) = b_e(:, :) end if if (present(x_e)) then ! Set Ax=b k = size(x_e, 1) l = size(x_e, 2) if (.not. allocated(obj%x_e)) then allocate (obj%x_e(k, l)) elseif (size(obj%x_e, 1) /= k .or. size(obj%x_e, 2) /= l) then deallocate (obj%x_e) allocate (obj%x_e(k, l)) end if obj%x_e(:, :) = x_e(:, :) end if end subroutine importLinearSolver !==================================================================================== !==================================================================================== subroutine prepareFixLinearSolver(obj, debug) class(LinearSolver_), intent(inout) :: obj logical, optional, intent(in) :: debug integer(int32), allocatable :: Index_I(:), Index_J(:), row_domain_id(:), column_Domain_ID(:) real(real64), allocatable:: val(:) integer(int32), allocatable :: array(:, :) integer(int32) :: i, m, n, rn, rd, cn, cd, same_n, count_reduc, j integer(int32) :: Index_I_max, Index_J_max, row_domain_id_max, column_Domain_ID_max ! It's too slow if (obj%ReadyForFix) return ! remove overlapped elements count_reduc = 0 if (present(debug)) then if (debug) then print *, "prepareFixLinearSolver >> [1] heap-sort started." end if end if ! first, heap sort n = size(obj%val) allocate (array(n, 4)) array(:, 1) = obj%Index_I array(:, 2) = obj%Index_J array(:, 3) = obj%row_domain_id array(:, 4) = obj%column_Domain_ID call heapsortArray(array, obj%val) obj%Index_I = array(:, 1) obj%Index_J = array(:, 2) obj%row_domain_id = array(:, 3) obj%column_Domain_ID = array(:, 4) if (present(debug)) then if (debug) then print *, "prepareFixLinearSolver >> [2] Remove overlaps" end if end if ! second, remove overlap do i = 1, n - 1 if (obj%index_i(i) == 0) cycle do j = i + 1, n if (obj%Index_I(i) == obj%Index_I(j) .and. obj%Index_J(i) == obj%Index_J(j)) then if (obj%row_domain_id(i) == obj%row_domain_id(j) .and. obj%column_domain_id(i) == obj%column_domain_id(j)) then obj%val(i) = obj%val(i) + obj%val(j) obj%Index_I(j) = 0 obj%row_domain_id(j) = 0 obj%Index_j(j) = 0 obj%column_domain_id(j) = 0 else exit end if else exit end if end do end do ! regacy ! do i=1,size(obj%Index_I)-1 ! if(obj%Index_I(i) == 0 ) cycle ! ! !$OMP parallel ! !$OMP do ! do j=i+1, size(obj%Index_I) ! if(obj%row_domain_id(j)==obj%row_domain_id(i) )then ! if(obj%column_domain_id(j)==obj%column_domain_id(i) )then ! if(obj%Index_I(j) == obj%Index_I(i) )then ! if(obj%Index_J(j) == obj%Index_J(i) )then ! obj%val(i) = obj%val(i) + obj%val(j) ! obj%Index_I(j) = 0 ! obj%row_domain_id(j) = 0 ! obj%Index_j(j) = 0 ! obj%column_domain_id(j) = 0 ! else ! cycle ! endif ! else ! cycle ! endif ! else ! cycle ! endif ! else ! cycle ! endif ! enddo ! !$OMP end do ! !$OMP end parallel ! enddo ! if (present(debug)) then if (debug) then print *, "prepareFixLinearSolver >> [3] Renew info" end if end if count_reduc = 0 do i = 1, size(obj%index_I) if (obj%index_I(i) /= 0) then count_reduc = count_reduc + 1 end if end do allocate (val(count_reduc)) allocate (Index_I(count_reduc)) allocate (Index_J(count_reduc)) allocate (row_domain_id(count_reduc)) allocate (column_Domain_ID(count_reduc)) n = 0 do i = 1, size(obj%Index_I) if (obj%Index_I(i) == 0) cycle n = n + 1 val(n) = obj%val(i) Index_I(n) = obj%Index_I(i) Index_J(n) = obj%Index_J(i) row_domain_id(n) = obj%row_domain_id(i) column_Domain_ID(n) = obj%column_Domain_ID(i) end do deallocate (obj%val) deallocate (obj%Index_I) deallocate (obj%Index_J) deallocate (obj%row_domain_id) deallocate (obj%column_Domain_ID) obj%val = val obj%Index_I = Index_I obj%Index_J = Index_J obj%row_domain_id = row_domain_id obj%column_Domain_ID = column_Domain_ID obj%ReadyForFix = .true. if (present(debug)) then if (debug) then print *, "prepareFixLinearSolver >> [ok] Done" end if end if end subroutine !==================================================================================== !==================================================================================== subroutine solveLinearSolver(obj, Solver, MPI, OpenCL, CUDAC, preconditioning, CRS) class(LinearSolver_), intent(inout) :: obj character(*), intent(in) :: Solver logical, optional, intent(in) :: MPI, OpenCL, CUDAC, preconditioning, CRS integer(int32), allocatable :: Index_I(:), Index_J(:), row_domain_id(:), column_Domain_ID(:) integer(int64), allocatable :: CRS_Index_I(:) real(real64), allocatable:: val(:) integer(int32) :: i, m, n, rn, rd, cn, cd, same_n, count_reduc, j print *, "[Warning] Please do not use this procedure" print *, "because this procedure has some bugs" print *, "Please use FEMSolverClass instead of this solver." ! if not allocated COO format if (.not. allocated(obj%val)) then if (allocated(obj%a) .and. allocated(obj%b)) then ! it only has obj%a and obj%b n = size(obj%b) obj%x = zeros(n) if (Solver == "BiCGSTAB") then call bicgstab1d(a=obj%a, b=obj%b, x=obj%x, n=n, itrmax=obj%itrmax, er=obj%er0) elseif (Solver == "GPBiCG") then call bicgstab1d(a=obj%a, b=obj%b, x=obj%x, n=n, itrmax=obj%itrmax, er=obj%er0) elseif (Solver == "GaussJordan") then call gauss_jordan_pv(obj%a, obj%x, obj%b, size(obj%b, 1)) else call gauss_jordan_pv(obj%a, obj%x, obj%b, size(obj%b, 1)) end if return end if end if if (.not. allocated(obj%a) .and. .not. allocated(obj%val)) then print *, "solveLinearSolver >> ERROR :: .not. allocated(obj%b) " stop end if if (.not. allocated(obj%b)) then print *, "solveLinearSolver >> ERROR :: .not. allocated(obj%b) " stop end if if (.not. allocated(obj%x)) then allocate (obj%x(size(obj%b))) obj%x(:) = 0.0d0 end if if (size(obj%x) /= size(obj%b)) then deallocate (obj%x) allocate (obj%x(size(obj%b))) obj%x(:) = 0.0d0 end if ! No MPI, No OpenCl and No CUDAC if (allocated(obj%a)) then if (allocated(obj%b)) then if (allocated(obj%x)) then ! run as non EBE-mode if (Solver == "GaussSeidel") then call gauss_seidel(obj%a, obj%b, obj%x, size(obj%a, 1), obj%itrmax, obj%er0) elseif (Solver == "GaussJordanPV" .or. Solver == "GaussJordan") then call gauss_jordan_pv(obj%a, obj%x, obj%b, size(obj%a, 1)) elseif (Solver == "BiCGSTAB") then if (present(CRS)) then if (CRS .eqv. .true.) then if (.not. allocated(obj%CRS_val)) then call obj%convertCOOtoCRS() end if call bicgstab_CRS(obj%val, obj%CRS_index_I, obj%CRS_index_J, obj%x, obj%b, obj%itrmax, obj%er0, obj%debug) elseif (allocated(obj%val)) then ! COO format if (allocated(obj%locked)) then call bicgstab_COO(obj%val, obj%index_I, obj%index_J, obj%x, obj%b, obj%itrmax, obj%er0, obj%debug, obj%locked) else call bicgstab_COO(obj%val, obj%index_I, obj%index_J, obj%x, obj%b, obj%itrmax, obj%er0, obj%debug) end if else call bicgstab1d(obj%a, obj%b, obj%x, size(obj%a, 1), obj%itrmax, obj%er0) end if else call bicgstab1d(obj%a, obj%b, obj%x, size(obj%a, 1), obj%itrmax, obj%er0) end if elseif (Solver == "GPBiCG") then if (present(preconditioning)) then if (preconditioning .eqv. .true.) then call preconditioned_GPBiCG(obj%a, obj%b, obj%x, size(obj%a, 1), obj%itrmax, obj%er0) else call GPBiCG(obj%a, obj%b, obj%x, size(obj%a, 1), obj%itrmax, obj%er0) end if else call GPBiCG(obj%a, obj%b, obj%x, size(obj%a, 1), obj%itrmax, obj%er0) end if else print *, "LinearSolver_ ERROR:: no such solver as :: ", Solver end if return end if end if else if (allocated(obj%NumberOfNode)) then ! May be overlapped! ! remove overlap print *, "solveLinearSolver >> preparing..." Index_I = obj%Index_I Index_J = obj%Index_J do i = 1, size(Index_I) m = obj%row_Domain_ID(i) if (m == 1) then n = 0 else n = sum(obj%NumberOfNode(1:m - 1))*obj%DOF end if Index_I(i) = Index_I(i) + n end do do i = 1, size(Index_J) m = obj%column_Domain_ID(i) if (m == 1) then n = 0 else n = sum(obj%NumberOfNode(1:m - 1))*obj%DOF end if Index_J(i) = Index_J(i) + n end do print *, "solveLinearSolver >> start!" if (present(CRS)) then ! create CRS_Index_I from Index_I@COO if (CRS) then call bicgstab_CRS(obj%val, CRS_index_I, index_J, obj%x, obj%b, obj%itrmax, obj%er0, obj%debug) return end if end if if (Solver == "BiCGSTAB") then if (allocated(obj%locked)) then call bicgstab_COO(obj%val, index_I, index_J, obj%x, obj%b, obj%itrmax, obj%er0, obj%debug, obj%locked) else call bicgstab_COO(obj%val, index_I, index_J, obj%x, obj%b, obj%itrmax, obj%er0, obj%debug) end if else call GaussJordan_COO(obj%val, index_I, index_J, obj%x, obj%b) end if else if (present(CRS)) then if (CRS) then call bicgstab_CRS(obj%val, CRS_index_I, index_J, obj%x, obj%b, obj%itrmax, obj%er0, obj%debug) return end if end if if (Solver == "BiCGSTAB") then if (allocated(obj%locked)) then call bicgstab_COO(obj%val, obj%index_I, obj%index_J, obj%x, obj%b, obj%itrmax, obj%er0, obj%debug & , obj%locked) else call bicgstab_COO(obj%val, obj%index_I, obj%index_J, obj%x, obj%b, obj%itrmax, obj%er0, obj%debug) end if else call GaussJordan_COO(obj%val, obj%index_I, obj%index_J, obj%x, obj%b) end if end if return end if print *, "LinearSolver_ ERROR:: EBE-mode is not implemented yet." stop end subroutine solveLinearSolver !==================================================================================== !==================================================================================== subroutine gauss_seidel(a, b, x, n, itrmax, er0) integer(int32), intent(in) :: n, itrmax real(real64), intent(in) :: a(n, n), b(n), er0 real(real64), intent(out) :: x(n) real(real64) s, er, rd(n), r(n) integer(int32) i, itr do i = 1, n if (a(i, i) == 0.0d0) stop 'a(i, i) == 0.0d0' rd(i) = 1.0d0/a(i, i) end do x(1:n) = 0.0d0 do itr = 1, itrmax do i = 1, n s = dot_product(a(i, 1:i - 1), x(1:i - 1)) s = s + dot_product(a(i, i + 1:n), x(i + 1:n)) x(i) = rd(i)*(b(i) - s) end do r(1:n) = b(1:n) - matmul(a, x) er = dot_product(r, r) if (er <= er0) then write (20, *) '# converged#' exit end if end do end subroutine gauss_seidel !=================================================================================== subroutine GaussJordan_COO(a0, index_i, index_j, x, b) integer(int32) :: n real(real64), intent(in) :: a0(:), b(:) integer(int32), intent(in) :: index_i(:), index_j(:) real(real64), allocatable, intent(inout) :: x(:) integer(int32) i, j, jj, k, m, nn, mm, id1, id2 real(real64), allocatable :: a(:, :), w(:) real(real32) ar, am, t print *, "Solver :: GaussJordan" ! CAUTION this code is not optimized. ! It just converts a COO-formatted sparse matrix ! to Dense matrix, and pass them to ! Gauss_jordan_pv_real64 n = size(b) a = zeros(n, n) x = zeros(n) print *, "DOF = ", n do i = 1, size(index_i) a(index_i(i), index_j(i)) = a0(i) end do nn = n x(:) = b(:) do k = 1, n m = k am = abs(a(k, k)) do i = k + 1, n if (abs(a(i, k)) > am) then am = abs(a(i, k)) m = i end if end do if (am == 0.0d0) stop ' A is singular ' if (k /= m) then w(k:n) = a(k, k:n) a(k, k:n) = a(m, k:n) a(m, k:n) = w(k:n) t = x(k) x(k) = x(m) x(m) = t end if ! gauss_jordan if (a(k, k) == 0.0d0) stop 'devide by zero3gauss_jordan_pv' ar = 1.0d0/a(k, k) a(k, k) = 1.0d0 a(k, k + 1:n) = ar*a(k, k + 1:n) x(k) = ar*x(k) do i = 1, n if (i /= k) then a(i, k + 1:n) = a(i, K + 1:n) - a(i, k)*a(k, k + 1:n) x(i) = x(i) - a(i, k)*x(k) a(i, k) = 0.0d0 end if end do end do return ! ! ! ! ! nn = size(b,1) ! n = size(b,1) ! a = a0 ! w = zeros(n) ! x = b ! ! do k = 1, n ! m = k ! ! ! am = A(k,k) ! am = 0.0d0 ! do j=1,size(index_i) ! if(index_i(j)==k .and. index_j(j)==k )then ! am = abs(a(j)) ! exit ! endif ! enddo ! ! ! ! !do i = k+1, n ! ! if (abs(a(i,k)) > am) then ! ! am = abs(a(i,k)) ! ! m = i ! ! endif ! !enddo ! ! do j=1,size(index_i) ! if(index_i(j)>=k+1 .and. k == index_j(j) )then ! if( abs(a(j)) > am )then ! am = abs(a(j) ) ! m = index_i(j) ! endif ! endif ! enddo ! ! ! ! if (am == 0.0d0) stop ' A is singular ' ! ! !if ( k /= m) then ! ! w(k:n) = a(k, k:n) ! ! a(k,k:n) = a(m, k:n) ! ! a(m, k:n) =w(k:n) ! ! t = x(k) ! ! x(k) = x(m) ! ! x(m) = t ! !endif ! if( k/=m)then ! !(1) w(k:n) = a(k, k:n) ! do j=1,size(index_i) ! id1 = index_i(j) ! if(id1==k)then ! id2 = index_j(j) ! if(k<=id2 .and. id2 <= n)then ! w(id2) = a(j) ! endif ! enddo ! enddo ! !(2) a(k,k:n) = a(m, k:n) ! do j=1,size(index_i) ! id1 = index_i(j) ! if(id1==k)then ! id2 = index_j(j) ! if(k<=id2 .and. id2 <= n)then ! w(id2) = a(j) ! endif ! enddo ! enddo ! ! ! ! endif ! ! ! ! gauss_jordan ! if (a(k, k) == 0.0d0) stop 'devide by zero3gauss_jordan_pv' ! ar = 1.0d0 / a(k,k) ! a(k,k) = 1.0d0 ! a(k,k+1:n) = ar * a(k, k+1:n) ! x(k) = ar * x(k) ! do i= 1, n ! if (i /= k) then ! a(i, k+1:n) = a(i, K+1:n) - a(i,k) * a(k, k+1:n) ! x(i) = x(i) - a(i,k) * x(k) ! a(i,k) = 0.0d0 ! endif ! enddo ! enddo end subroutine !=========================================================================== !=================================================================================== subroutine gauss_jordan_pv_real64(a0, x, b, n) integer(int32), intent(in) :: n real(real64), intent(in) :: a0(n, n), b(n) real(real64), intent(inout) :: x(n) integer(int32) i, j, k, m, nn, mm real(real64) ar, am, t, a(n, n), w(n) nn = size(a0, 1) a(:, :) = a0(:, :) x(:) = b(:) do k = 1, n m = k am = abs(a(k, k)) do i = k + 1, n if (abs(a(i, k)) > am) then am = abs(a(i, k)) m = i end if end do if (am == 0.0d0) stop ' A is singular ' if (k /= m) then w(k:n) = a(k, k:n) a(k, k:n) = a(m, k:n) a(m, k:n) = w(k:n) t = x(k) x(k) = x(m) x(m) = t end if ! gauss_jordan if (a(k, k) == 0.0d0) stop 'devide by zero3gauss_jordan_pv' ar = 1.0d0/a(k, k) a(k, k) = 1.0d0 a(k, k + 1:n) = ar*a(k, k + 1:n) x(k) = ar*x(k) do i = 1, n if (i /= k) then a(i, k + 1:n) = a(i, K + 1:n) - a(i, k)*a(k, k + 1:n) x(i) = x(i) - a(i, k)*x(k) a(i, k) = 0.0d0 end if end do end do end subroutine !=========================================================================== !=================================================================================== subroutine gauss_jordan_pv_real32(a0, x, b, n) integer(int32), intent(in) :: n real(real32), intent(in) :: a0(n, n), b(n) real(real32), intent(out) :: x(n) integer(int32) i, j, k, m, nn, mm real(real32) ar, am, t, a(n, n), w(n) nn = size(a0, 1) a(:, :) = a0(:, :) x(:) = b(:) do k = 1, n m = k am = abs(a(k, k)) do i = k + 1, n if (abs(a(i, k)) > am) then am = abs(a(i, k)) m = i end if end do if (am == 0.0d0) stop ' A is singular ' if (k /= m) then w(k:n) = a(k, k:n) a(k, k:n) = a(m, k:n) a(m, k:n) = w(k:n) t = x(k) x(k) = x(m) x(m) = t end if ! gauss_jordan if (a(k, k) == 0.0d0) stop 'devide by zero3gauss_jordan_pv' ar = 1.0d0/a(k, k) a(k, k) = 1.0d0 a(k, k + 1:n) = ar*a(k, k + 1:n) x(k) = ar*x(k) do i = 1, n if (i /= k) then a(i, k + 1:n) = a(i, K + 1:n) - a(i, k)*a(k, k + 1:n) x(i) = x(i) - a(i, k)*x(k) a(i, k) = 0.0d0 end if end do end do end subroutine !=========================================================================== !=================================================================================== subroutine gauss_jordan_pv_complex64(a0, x, b, n) integer(int32), intent(in) :: n complex(real64), intent(in) :: a0(n, n), b(n) complex(real64), intent(out) :: x(n) integer(int32) i, j, k, m, nn, mm complex(real64) ar, am, t, a(n, n), w(n) nn = size(a0, 1) a(:, :) = a0(:, :) x(:) = b(:) do k = 1, n m = k am = abs(a(k, k)) do i = k + 1, n if (abs(a(i, k)) > abs(am)) then am = abs(a(i, k)) m = i end if end do if (am == 0.0d0) stop ' A is singular ' if (k /= m) then w(k:n) = a(k, k:n) a(k, k:n) = a(m, k:n) a(m, k:n) = w(k:n) t = x(k) x(k) = x(m) x(m) = t end if ! gauss_jordan if (a(k, k) == 0.0d0) stop 'devide by zero3gauss_jordan_pv' ar = 1.0d0/a(k, k) a(k, k) = 1.0d0 a(k, k + 1:n) = ar*a(k, k + 1:n) x(k) = ar*x(k) do i = 1, n if (i /= k) then a(i, k + 1:n) = a(i, K + 1:n) - a(i, k)*a(k, k + 1:n) x(i) = x(i) - a(i, k)*x(k) a(i, k) = 0.0d0 end if end do end do end subroutine !=========================================================================== subroutine bicgstab_diffusion(a, b, x, n, itrmax, er, DBC, DBCVal) integer(int32), intent(in) :: n, itrmax, DBC(:, :) real(real64), intent(in) :: a(n, n), b(n), er, DBCVal(:, :) real(real64), intent(inout) :: x(n) integer(int32) itr, i, j real(real64) alp, bet, c1, c2, c3, ev, vv, rr, er0, init_rr real(real64) r(n), r0(n), p(n), y(n), e(n), v(n) er0 = 1.00e-8 r(:) = b - matmul(a, x) ! if DBC => reset residual do i = 1, size(DBC, 1) do j = 1, size(DBC, 2) if (DBC(i, j) < 1) then cycle else r(DBC(i, j)) = 0.0d0 x(DBC(i, j)) = DBCVal(i, j) end if end do end do ! c1 = dot_product(r, r) init_rr = c1 if (c1 == 0.0d0) then print *, "Caution :: Initial residual is zero" return end if p(:) = r(:) r0(:) = r(:) do itr = 1, itrmax !c1 = dot_product(r0,r) y(:) = matmul(a, p) c2 = dot_product(r0, y) alp = c1/c2 e(:) = r(:) - alp*y(:) v(:) = matmul(a, e) ev = dot_product(e, v) vv = dot_product(v, v) if (vv == 0.0d0) stop "Bicgstab devide by zero" c3 = ev/vv x(:) = x(:) + alp*p(:) + c3*e(:) r(:) = e(:) - c3*v(:) ! if DBC => reset residual do i = 1, size(DBC, 1) do j = 1, size(DBC, 2) if (DBC(i, j) < 1) then cycle else r(DBC(i, j)) = 0.0d0 x(DBC(i, j)) = DBCVal(i, j) end if end do end do rr = dot_product(r, r) ! write(*,*) 'itr, er =', itr,rr if (rr/init_rr < er0) then print *, "[ok] :: BICGSTAB is converged in ", i, " steps." exit end if c1 = dot_product(r0, r) bet = c1/(c2*c3) if ((c2*c3) == 0.0d0) stop "Bicgstab devide by zero" p(:) = r(:) + bet*(p(:) - c3*y(:)) if (itrmax == itr) then print *, "ERROR :: BICGSTAB did not converge" return end if end do end subroutine bicgstab_diffusion !=============================================================== subroutine bicgstab_CRS(a, ptr_i, index_j, x, b, itrmax, er, debug) integer(int32), intent(inout) :: index_j(:), itrmax integer(int64), intent(inout) :: ptr_i(:) real(real64), intent(inout) :: a(:), b(:), er real(real64), intent(inout) :: x(:) logical, optional, intent(in) :: debug logical :: speak = .false. integer(int32) itr, i, j, n real(real64) alp, bet, c1, c2, c3, ev, vv, rr, er0, init_rr real(real64), allocatable:: r(:), r0(:), p(:), y(:), e(:), v(:), pa(:), ax(:) er0 = er if (present(debug)) then speak = debug end if n = size(b) if (speak) print *, "BiCGSTAB STARTED >> DOF:", n allocate (r(n), r0(n), p(n), y(n), e(n), v(n)) r(:) = b(:) if (speak) print *, "BiCGSTAB >> [1] initialize" ax = crs_matvec(CRS_value=a, CRS_col=index_j, & CRS_row_ptr=ptr_i, old_vector=x) r = b - ax if (speak) print *, "BiCGSTAB >> [2] dp1" c1 = dot_product(r, r) init_rr = c1 if (speak) print *, "BiCGSTAB >> |r|^2 = ", init_rr if (c1 < er0) return p(:) = r(:) r0(:) = r(:) do itr = 1, itrmax if (speak) print *, "BiCGSTAB >> ["//str(itr)//"] initialize" c1 = dot_product(r0, r) y = crs_matvec(CRS_value=a, CRS_col=index_j, & CRS_row_ptr=ptr_i, old_vector=p) c2 = dot_product(r0, y) alp = c1/c2 e(:) = r(:) - alp*y(:) v = crs_matvec(CRS_value=a, CRS_col=index_j, & CRS_row_ptr=ptr_i, old_vector=e) if (speak) print *, "BiCGSTAB >> ["//str(itr)//"] half" ev = dot_product(e, v) vv = dot_product(v, v) if (vv == 0.0d0) stop "Bicgstab devide by zero" c3 = ev/vv x(:) = x(:) + alp*p(:) + c3*e(:) r(:) = e(:) - c3*v(:) rr = dot_product(r, r) if (speak) then print *, rr/init_rr end if ! write(*,*) 'itr, er =', itr,rr if (rr/init_rr < er0) exit c1 = dot_product(r0, r) bet = c1/(c2*c3) if ((c2*c3) == 0.0d0) stop "Bicgstab devide by zero" p(:) = r(:) + bet*(p(:) - c3*y(:)) end do end subroutine !=============================================================== !=============================================================== subroutine bicgstab_COO(a, index_i, index_j, x, b, itrmax, er, debug, locked) integer(int32), intent(inout) :: index_i(:), index_j(:), itrmax real(real64), intent(inout) :: a(:), b(:), er real(real64), intent(inout) :: x(:) logical, optional, intent(in) :: debug, locked(:) logical :: speak = .false. integer(int32) itr, i, j, n real(real64) alp, bet, c1, c2, c3, ev, vv, rr, er0, init_rr real(real64), allocatable:: r(:), r0(:), p(:), y(:), e(:), v(:) print *, "[ok] BiCGSTAB for COO started." if (present(debug)) then speak = debug end if if (speak) print *, "BiCGSTAB STARTED >> DOF:", n n = size(b) allocate (r(n), r0(n), p(n), y(n), e(n), v(n)) er0 = er r(:) = b(:) if (speak) print *, "BiCGSTAB >> [1] initialize" do i = 1, size(a) if (index_i(i) <= 0) cycle r(index_i(i)) = r(index_i(i)) - a(i)*x(index_j(i)) end do if (speak) print *, dot_product(r, r), size(a) ! >> if fixed >> r=0, x=b if (present(locked)) then !!$OMP parallel do private(i) do i = 1, size(locked) if (locked(i)) then r(i) = 0.0d0 x(i) = b(i) end if end do !!$OMP end parallel do end if !r(:) = b - matmul(a,x) c1 = dot_product(r, r) if (speak) print *, "BiCGSTAB >> [2] dp1", c1 init_rr = c1 if (c1 < er0) return p(:) = r(:) r0(:) = r(:) do itr = 1, itrmax if (speak) print *, "BiCGSTAB >> ["//str(itr)//"] initialize" c1 = dot_product(r0, r) !y(:) = matmul(a,p) y(:) = 0.0d0 !!$OMP parallel do reduction(+:y) private(i) do i = 1, size(a) if (index_i(i) <= 0) then ! do nothing else y(index_i(i)) = y(index_i(i)) + a(i)*p(index_j(i)) end if end do !!$OMP end parallel do c2 = dot_product(r0, y) alp = c1/c2 e(:) = r(:) - alp*y(:) !v(:) = matmul(a,e) if (speak) print *, "BiCGSTAB >> ["//str(itr)//"] half" v(:) = 0.0d0 !!$OMP parallel do reduction(+:v) private(i) do i = 1, size(a) if (index_i(i) <= 0) then ! do nothing else v(index_i(i)) = v(index_i(i)) + a(i)*e(index_j(i)) end if end do !!$OMP end parallel do ev = dot_product(e, v) vv = dot_product(v, v) if (vv == 0.0d0) stop "Bicgstab devide by zero" c3 = ev/vv x(:) = x(:) + alp*p(:) + c3*e(:) r(:) = e(:) - c3*v(:) ! >> if fixed >> r=0, x=b if (present(locked)) then !!$OMP parallel do private(i) do i = 1, size(locked) if (locked(i)) then r(i) = 0.0d0 x(i) = b(i) end if end do !! $OMP end parallel do end if rr = dot_product(r, r) if (speak) print *, 'itr, |er|/|er0| =', itr, rr/init_rr if (rr/init_rr < er0) exit c1 = dot_product(r0, r) bet = c1/(c2*c3) if ((c2*c3) == 0.0d0) stop "Bicgstab devide by zero" p(:) = r(:) + bet*(p(:) - c3*y(:)) end do print *, "[ok] BiCGSTAB_COO >> error norm", rr/init_rr end subroutine !=============================================================== !=============================================================== subroutine bicgstab_real64(a, b, x, n, itrmax, er) integer(int32), intent(in) :: n, itrmax real(real64), intent(in) :: a(n, n), b(n), er real(real64), intent(inout) :: x(n) integer(int32) itr real(real64) alp, bet, c1, c2, c3, ev, vv, rr, er0, init_rr real(real64) r(n), r0(n), p(n), y(n), e(n), v(n) er0 = er r(:) = b - matmul(a, x) c1 = dot_product(r, r) init_rr = c1 if (c1 < er0) return p(:) = r(:) r0(:) = r(:) do itr = 1, itrmax c1 = dot_product(r0, r) y(:) = matmul(a, p) c2 = dot_product(r0, y) alp = c1/c2 e(:) = r(:) - alp*y(:) v(:) = matmul(a, e) ev = dot_product(e, v) vv = dot_product(v, v) if (vv == 0.0d0) stop "Bicgstab devide by zero" c3 = ev/vv x(:) = x(:) + alp*p(:) + c3*e(:) r(:) = e(:) - c3*v(:) rr = dot_product(r, r) ! write(*,*) 'itr, er =', itr,rr if (rr/init_rr < er0) exit c1 = dot_product(r0, r) bet = c1/(c2*c3) if ((c2*c3) == 0.0d0) stop "Bicgstab devide by zero" p(:) = r(:) + bet*(p(:) - c3*y(:)) end do end subroutine !=============================================================== !=============================================================== subroutine bicgstab_real32(a, b, x, n, itrmax, er) integer(int32), intent(in) :: n, itrmax real(real32), intent(in) :: a(n, n), b(n), er real(real32), intent(inout) :: x(n) integer(int32) itr real(real32) alp, bet, c1, c2, c3, ev, vv, rr, er0, init_rr real(real32) r(n), r0(n), p(n), y(n), e(n), v(n) er0 = er r(:) = b - matmul(a, x) c1 = dot_product(r, r) init_rr = c1 if (c1 < er0) return p(:) = r(:) r0(:) = r(:) do itr = 1, itrmax c1 = dot_product(r0, r) y(:) = matmul(a, p) c2 = dot_product(r0, y) alp = c1/c2 e(:) = r(:) - alp*y(:) v(:) = matmul(a, e) ev = dot_product(e, v) vv = dot_product(v, v) if (vv == 0.0d0) stop "Bicgstab devide by zero" c3 = ev/vv x(:) = x(:) + alp*p(:) + c3*e(:) r(:) = e(:) - c3*v(:) rr = dot_product(r, r) ! write(*,*) 'itr, er =', itr,rr if (rr/init_rr < er0) exit c1 = dot_product(r0, r) bet = c1/(c2*c3) if ((c2*c3) == 0.0d0) stop "Bicgstab devide by zero" p(:) = r(:) + bet*(p(:) - c3*y(:)) end do end subroutine !=============================================================== !=============================================================== subroutine bicgstab_complex64(a, b, x, n, itrmax, er) integer(int32), intent(in) :: n, itrmax complex(real64), intent(in) :: a(n, n), b(n), er complex(real64), intent(inout) :: x(n) integer(int32) itr complex(real64) alp, bet, c1, c2, c3, ev, vv, rr, er0, init_rr complex(real64) r(n), r0(n), p(n), y(n), e(n), v(n) er0 = er r(:) = b - matmul(a, x) c1 = dot_product(r, r) init_rr = c1 if (abs(c1) < abs(er0)) return p(:) = r(:) r0(:) = r(:) do itr = 1, itrmax c1 = dot_product(r0, r) y(:) = matmul(a, p) c2 = dot_product(r0, y) alp = c1/c2 e(:) = r(:) - alp*y(:) v(:) = matmul(a, e) ev = dot_product(e, v) vv = dot_product(v, v) if (vv == 0.0d0) stop "Bicgstab devide by zero" c3 = ev/vv x(:) = x(:) + alp*p(:) + c3*e(:) r(:) = e(:) - c3*v(:) rr = dot_product(r, r) ! write(*,*) 'itr, er =', itr,rr if (abs(rr/init_rr) < abs(er0)) exit c1 = dot_product(r0, r) bet = c1/(c2*c3) if ((c2*c3) == 0.0d0) stop "Bicgstab devide by zero" p(:) = r(:) + bet*(p(:) - c3*y(:)) end do end subroutine !=============================================================== !=============================================================== subroutine bicgstab1d(a, b, x, n, itrmax, er) integer(int32), intent(in) :: n, itrmax real(real64), intent(in) :: a(n, n), b(n), er real(real64), intent(inout) :: x(n) integer(int32) itr real(real64) alp, bet, c1, c2, c3, ev, vv, rr, er0, init_rr real(real64) r(n), r0(n), p(n), y(n), e(n), v(n) er0 = er r(:) = b - matmul(a, x) c1 = dot_product(r, r) init_rr = c1 if (c1 < er0) return p(:) = r(:) r0(:) = r(:) do itr = 1, itrmax c1 = dot_product(r0, r) y(:) = matmul(a, p) c2 = dot_product(r0, y) alp = c1/c2 e(:) = r(:) - alp*y(:) v(:) = matmul(a, e) ev = dot_product(e, v) vv = dot_product(v, v) if (vv == 0.0d0) stop "Bicgstab devide by zero" c3 = ev/vv x(:) = x(:) + alp*p(:) + c3*e(:) r(:) = e(:) - c3*v(:) rr = dot_product(r, r) ! write(*,*) 'itr, er =', itr,rr if (rr/init_rr < er0) exit c1 = dot_product(r0, r) bet = c1/(c2*c3) if ((c2*c3) == 0.0d0) stop "Bicgstab devide by zero" p(:) = r(:) + bet*(p(:) - c3*y(:)) end do end subroutine bicgstab1d !=============================================================== subroutine bicgstab_nr(a, b, x, n, itrmax, er, u_nod_x, u_nod_y) integer(int32), intent(in) :: n, itrmax, u_nod_x(:), u_nod_y(:) real(real64), intent(in) :: a(n, n), b(n), er real(real64), intent(inout) :: x(n) integer(int32) itr real(real64) alp, bet, c1, c2, c3, ev, vv, rr, er0, init_rr real(real64) r(n), r0(n), p(n), y(n), e(n), v(n) er0 = 1.00e-14 r(:) = b - matmul(a, x) call modify_residual(r, u_nod_x, u_nod_y) c1 = dot_product(r, r) init_rr = c1 if (c1 < er0) return p(:) = r(:) r0(:) = r(:) do itr = 1, itrmax y(:) = matmul(a, p) c2 = dot_product(r0, y) alp = c1/c2 e(:) = r(:) - alp*y(:) v(:) = matmul(a, e) ev = dot_product(e, v) vv = dot_product(v, v) if (vv == 0.0d0) stop "Bicgstab devide by zero" c3 = ev/vv x(:) = x(:) + alp*p(:) + c3*e(:) r(:) = e(:) - c3*v(:) call modify_residual(r, u_nod_x, u_nod_y) rr = dot_product(r, r) ! write(*,*) 'itr, er =', itr,rr if (rr/init_rr < er0) exit c1 = dot_product(r0, r) bet = c1/(c2*c3) if ((c2*c3) == 0.0d0) stop "Bicgstab devide by zero" p(:) = r(:) + bet*(p(:) - c3*y(:)) end do end subroutine bicgstab_nr !==================================================================================== subroutine bicgstab_nr1(a, b, x, n, itrmax, er, u_nod_x, u_nod_y, u_nod_dis_x, u_nod_dis_y) integer(int32), intent(in) :: n, itrmax, u_nod_x(:), u_nod_y(:) real(real64), intent(in) :: a(n, n), b(n), er, u_nod_dis_x(:), u_nod_dis_y(:) real(real64), intent(inout) :: x(n) integer(int32) itr real(real64) alp, bet, c1, c2, c3, ev, vv, rr, er0, init_rr real(real64) r(n), r0(n), p(n), y(n), e(n), v(n) er0 = 1.00e-14 r(:) = b - matmul(a, x) call modify_residual_1(r, x, u_nod_x, u_nod_y, u_nod_dis_x, u_nod_dis_y) c1 = dot_product(r, r) init_rr = c1 if (c1 < er0) return p(:) = r(:) r0(:) = r(:) do itr = 1, itrmax y(:) = matmul(a, p) c2 = dot_product(r0, y) alp = c1/c2 e(:) = r(:) - alp*y(:) v(:) = matmul(a, e) ev = dot_product(e, v) vv = dot_product(v, v) if (vv == 0.0d0) stop "Bicgstab devide by zero" c3 = ev/vv x(:) = x(:) + alp*p(:) + c3*e(:) r(:) = e(:) - c3*v(:) call modify_residual_1(r, x, u_nod_x, u_nod_y, u_nod_dis_x, u_nod_dis_y) rr = dot_product(r, r) ! write(*,*) 'itr, er =', itr,rr if (rr/init_rr < er0) exit c1 = dot_product(r0, r) bet = c1/(c2*c3) if ((c2*c3) == 0.0d0) stop "Bicgstab devide by zero" p(:) = r(:) + bet*(p(:) - c3*y(:)) end do end subroutine bicgstab_nr1 !==================================================================================== subroutine bicgstab_dirichlet(a, b, x, n, itrmax, er, DBoundNodID, DBoundVal, SetBC) integer(int32), intent(in) :: n, itrmax, DBoundNodID(:, :), SetBC real(real64), intent(in) :: a(n, n), b(n), er, DBoundVal(:, :) real(real64), intent(inout) :: x(n) integer(int32) itr real(real64) alp, bet, c1, c2, c3, ev, vv, rr, er0, init_rr real(real64) r(n), r0(n), p(n), y(n), e(n), v(n) er0 = 1.00e-14 r(:) = b - matmul(a, x) call modify_residual_dirichlet(r, x, DBoundNodID, DBoundVal, SetBC) c1 = dot_product(r, r) init_rr = c1 if (c1 < er0) return p(:) = r(:) r0(:) = r(:) do itr = 1, itrmax y(:) = matmul(a, p) c2 = dot_product(r0, y) alp = c1/c2 e(:) = r(:) - alp*y(:) v(:) = matmul(a, e) ev = dot_product(e, v) vv = dot_product(v, v) if (vv == 0.0d0) stop "Bicgstab devide by zero" c3 = ev/vv x(:) = x(:) + alp*p(:) + c3*e(:) r(:) = e(:) - c3*v(:) call modify_residual_dirichlet(r, x, DBoundNodID, DBoundVal, SetBC) rr = dot_product(r, r) ! write(*,*) 'itr, er =', itr,rr if (rr/init_rr < er0) exit c1 = dot_product(r0, r) bet = c1/(c2*c3) if ((c2*c3) == 0.0d0) stop "Bicgstab devide by zero" p(:) = r(:) + bet*(p(:) - c3*y(:)) end do end subroutine bicgstab_dirichlet !==================================================================================== subroutine modify_residual_1(r, x, u_nod_x, u_nod_y, u_nod_dis_x, u_nod_dis_y) integer(int32), intent(in)::u_nod_x(:), u_nod_y(:) real(real64), intent(in) :: u_nod_dis_x(:), u_nod_dis_y(:) real(real64), intent(inout)::r(:), x(:) integer(int32) i do i = 1, size(u_nod_x) r(2*u_nod_x(i) - 1) = 0.0d0 x(2*u_nod_x(i) - 1) = u_nod_dis_x(i) end do do i = 1, size(u_nod_y) r(2*u_nod_y(i)) = 0.0d0 x(2*u_nod_y(i)) = u_nod_dis_y(i) end do end subroutine modify_residual_1 !==================================================================================== subroutine modify_residual(r, u_nod_x, u_nod_y) integer(int32), intent(in)::u_nod_x(:), u_nod_y(:) real(real64), intent(inout)::r(:) integer(int32) i do i = 1, size(u_nod_x) r(2*u_nod_x(i) - 1) = 0.0d0 end do do i = 1, size(u_nod_y) r(2*u_nod_y(i)) = 0.0d0 end do end subroutine modify_residual !==================================================================================== subroutine modify_residual_dirichlet(r, x, DBoundNodID, DBoundVal, SetBC) integer(int32), intent(in)::DBoundNodID(:, :), SetBC real(real64), intent(in)::DBoundVal(:, :) real(real64), intent(inout)::r(:), x(:) integer(int32) :: i, j, k, dim_num, dbc_num real(real64) :: val if (SetBC == 1) then dim_num = size(DBoundNodID, 2) dbc_num = size(DBoundNodID, 1) do i = 1, dim_num do j = 1, dbc_num k = DBoundNodID(j, i) val = DBoundVal(j, i) if (k < 1) then cycle elseif (k >= 1) then x(dim_num*(k - 1) + i) = val*dble(SetBC) r(dim_num*(k - 1) + i) = 0.0d0 else cycle end if end do end do else dim_num = size(DBoundNodID, 2) dbc_num = size(DBoundNodID, 1) do i = 1, dim_num do j = 1, dbc_num k = DBoundNodID(j, i) val = DBoundVal(j, i) if (k < 1) then cycle elseif (k >= 1) then r(dim_num*(k - 1) + i) = 0.0d0 x(dim_num*(k - 1) + i) = 0.0d0 else cycle end if end do end do end if end subroutine !==================================================================================== !=========================================================================== subroutine GPBiCG_real64(a, b, x, n, itrmax, er) integer(int32), intent(in) :: n integer(int32), optional, intent(in) :: itrmax real(real64), intent(in) :: a(n, n), b(n) real(real64), optional, intent(in)::er real(real64), intent(inout) :: x(n) integer(int32) itr, itrmax_ real(real64) alp, c1, rr, er0, init_rr, beta real(real64) gzi, nu, val1, val2, r0rk, eps real(real64) r(n), r0(n), p(n), y(n), ap(n), q(n) real(real64) u(n), w(n), t(n), t_(n), z(n) eps = 1.0e-18 er0 = input(default=eps, option=er) r(:) = b - matmul(a, x) c1 = dot_product(r, r) init_rr = c1 if (c1 < er0) return beta = 0.0d0 p(:) = 0.0d0 r0(:) = r(:) w(:) = 0.0d0 u(:) = 0.0d0 t(:) = 0.0d0 t_(:) = 0.0d0 z(:) = 0.0d0 itrmax_ = input(default=1000, option=itrmax) do itr = 1, itrmax_ p(:) = r(:) + beta*(p(:) - u(:)) ! triple checked. ap(:) = matmul(a, p) ! triple checked. ap=v, alp = dot_product(r0, r)/dot_product(r0, ap)! triple checked. y(:) = t(:) - r(:) - alp*w(:) + alp*ap(:) !triple checked. t(:) = r(:) - alp*ap(:) ! triple checked. q(:) = matmul(a, t) ! triple checked. s=q=at=c r0rk = dot_product(r0, r) if (itr == 1) then gzi = dot_product(q, t)/dot_product(q, q)! double checked. nu = 0.0d0 ! double checked. else val1 = dot_product(y, y)*dot_product(q, t) - dot_product(y, t)*dot_product(q, y) val2 = dot_product(q, q)*dot_product(y, y) - dot_product(y, q)*dot_product(q, y) if (val2 == 0.0d0) then !print *, itr !print *, r(:), alp*ap(:) stop "GPBiCG devide by zero" end if gzi = val1/val2 val1 = dot_product(q, q)*dot_product(y, t) - dot_product(y, q)*dot_product(q, t) nu = val1/val2 !triple checked. end if u(:) = gzi*ap + nu*(t_(:) - r(:) + beta*u(:)) !double checked. z(:) = gzi*r(:) + nu*z(:) - alp*u(:)!double checked. x(:) = x(:) + alp*p(:) + z(:) !double checked. r(:) = t(:) - nu*y(:) - gzi*q(:)!double checked. t_(:) = t(:) rr = dot_product(r, r) !print *, 'itr, er =', itr,sqrt(rr),sqrt(rr)/sqrt(init_rr) !print *, "ans = ",x(:) if (sqrt(rr)/sqrt(init_rr) < er0) exit r0rk = dot_product(r0, r) beta = alp/gzi*dot_product(r0, r)/r0rk w(:) = q(:) + beta*ap(:) if (itr == itrmax_) then print *, "ERROR :: GPBiCG did not converge." end if end do end subroutine !=============================================================== !=========================================================================== subroutine GPBiCG_real32(a, b, x, n, itrmax, er) integer(int32), intent(in) :: n integer(int32), optional, intent(in) :: itrmax real(real32), intent(in) :: a(n, n), b(n) real(real32), optional, intent(in)::er real(real32), intent(inout) :: x(n) integer(int32) itr, itrmax_ real(real32) alp, c1, rr, er0, init_rr, beta real(real32) gzi, nu, val1, val2, r0rk, eps real(real32) r(n), r0(n), p(n), y(n), ap(n), q(n) real(real32) u(n), w(n), t(n), t_(n), z(n) eps = 1.0e-18 er0 = input(default=eps, option=er) r(:) = b - matmul(a, x) c1 = dot_product(r, r) init_rr = c1 if (c1 < er0) return beta = 0.0d0 p(:) = 0.0d0 r0(:) = r(:) w(:) = 0.0d0 u(:) = 0.0d0 t(:) = 0.0d0 t_(:) = 0.0d0 z(:) = 0.0d0 itrmax_ = input(default=1000, option=itrmax) do itr = 1, itrmax_ p(:) = r(:) + beta*(p(:) - u(:)) ! triple checked. ap(:) = matmul(a, p) ! triple checked. ap=v, alp = dot_product(r0, r)/dot_product(r0, ap)! triple checked. y(:) = t(:) - r(:) - alp*w(:) + alp*ap(:) !triple checked. t(:) = r(:) - alp*ap(:) ! triple checked. q(:) = matmul(a, t) ! triple checked. s=q=at=c r0rk = dot_product(r0, r) if (itr == 1) then gzi = dot_product(q, t)/dot_product(q, q)! double checked. nu = 0.0d0 ! double checked. else val1 = dot_product(y, y)*dot_product(q, t) - dot_product(y, t)*dot_product(q, y) val2 = dot_product(q, q)*dot_product(y, y) - dot_product(y, q)*dot_product(q, y) if (val2 == 0.0d0) then !print *, itr !print *, r(:), alp*ap(:) stop "GPBiCG devide by zero" end if gzi = val1/val2 val1 = dot_product(q, q)*dot_product(y, t) - dot_product(y, q)*dot_product(q, t) nu = val1/val2 !triple checked. end if u(:) = gzi*ap + nu*(t_(:) - r(:) + beta*u(:)) !double checked. z(:) = gzi*r(:) + nu*z(:) - alp*u(:)!double checked. x(:) = x(:) + alp*p(:) + z(:) !double checked. r(:) = t(:) - nu*y(:) - gzi*q(:)!double checked. t_(:) = t(:) rr = dot_product(r, r) !print *, 'itr, er =', itr,sqrt(rr),sqrt(rr)/sqrt(init_rr) !print *, "ans = ",x(:) if (sqrt(rr)/sqrt(init_rr) < er0) exit r0rk = dot_product(r0, r) beta = alp/gzi*dot_product(r0, r)/r0rk w(:) = q(:) + beta*ap(:) if (itr == itrmax_) then print *, "ERROR :: GPBiCG did not converge." end if end do end subroutine !=============================================================== !=========================================================================== subroutine GPBiCG_complex64(a, b, x, n, itrmax, er) integer(int32), intent(in) :: n integer(int32), optional, intent(in) :: itrmax complex(real64), intent(in) :: a(n, n), b(n) complex(real64), optional, intent(in)::er complex(real64), intent(inout) :: x(n) integer(int32) itr, itrmax_ complex(real64) alp, c1, rr, er0, init_rr, beta complex(real64) gzi, nu, val1, val2, r0rk, eps complex(real64) r(n), r0(n), p(n), y(n), ap(n), q(n) complex(real64) u(n), w(n), t(n), t_(n), z(n) eps = 1.0e-18 er0 = input(default=eps, option=er) r(:) = b - matmul(a, x) c1 = dot_product(r, r) init_rr = c1 if (abs(c1) < abs(er0)) return beta = 0.0d0 p(:) = 0.0d0 r0(:) = r(:) w(:) = 0.0d0 u(:) = 0.0d0 t(:) = 0.0d0 t_(:) = 0.0d0 z(:) = 0.0d0 itrmax_ = input(default=1000, option=itrmax) do itr = 1, itrmax_ p(:) = r(:) + beta*(p(:) - u(:)) ! triple checked. ap(:) = matmul(a, p) ! triple checked. ap=v, alp = dot_product(r0, r)/dot_product(r0, ap)! triple checked. y(:) = t(:) - r(:) - alp*w(:) + alp*ap(:) !triple checked. t(:) = r(:) - alp*ap(:) ! triple checked. q(:) = matmul(a, t) ! triple checked. s=q=at=c r0rk = dot_product(r0, r) if (itr == 1) then gzi = dot_product(q, t)/dot_product(q, q)! double checked. nu = 0.0d0 ! double checked. else val1 = dot_product(y, y)*dot_product(q, t) - dot_product(y, t)*dot_product(q, y) val2 = dot_product(q, q)*dot_product(y, y) - dot_product(y, q)*dot_product(q, y) if (val2 == 0.0d0) then !print *, itr !print *, r(:), alp*ap(:) stop "GPBiCG devide by zero" end if gzi = val1/val2 val1 = dot_product(q, q)*dot_product(y, t) - dot_product(y, q)*dot_product(q, t) nu = val1/val2 !triple checked. end if u(:) = gzi*ap + nu*(t_(:) - r(:) + beta*u(:)) !double checked. z(:) = gzi*r(:) + nu*z(:) - alp*u(:)!double checked. x(:) = x(:) + alp*p(:) + z(:) !double checked. r(:) = t(:) - nu*y(:) - gzi*q(:)!double checked. t_(:) = t(:) rr = dot_product(r, r) !print *, 'itr, er =', itr,sqrt(rr),sqrt(rr)/sqrt(init_rr) !print *, "ans = ",x(:) if (abs(sqrt(rr)/sqrt(init_rr)) < abs(er0)) exit r0rk = dot_product(r0, r) beta = alp/gzi*dot_product(r0, r)/r0rk w(:) = q(:) + beta*ap(:) if (itr == itrmax_) then print *, "ERROR :: GPBiCG did not converge." end if end do end subroutine !=============================================================== !=========================================================================== subroutine preconditioned_GPBiCG(a, b, x, n, itrmax, er) integer(int32), intent(in) :: n integer(int32), optional, intent(in) :: itrmax real(real64), intent(in) :: a(1:n, 1:n), b(1:n) real(real64), optional, intent(in)::er real(real64), intent(inout) :: x(1:n) real(real64) :: L(1:n, 1:n), d(1:n) ! presented by Moe Thuthu et al., 2009, algorithm #3 integer(int32) itr, itrmax_, i, j, k real(real64) alp, c1, rr, er0, init_rr, beta, lld, ld real(real64) gzi, nu, val1, val2, r0rk, eps real(real64) r(n), rk(n), r_(n), r_k(n), r0(n), p(n), y(n), ap(n), q(n) real(real64) u(n), w(n), t(n), t_(n), z(n) print *, "<<< Under implementation >>>" ! Incomplete Cholosky decomposition. ! http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?%A5%B3%A5%EC%A5%B9%A5%AD%A1%BC%CA%AC%B2%F2 ! >>>>>>>>>> L(:, :) = 0.0d0 d(:) = 0.0d0 d(1) = a(1, 1) L(1, 1) = 1.0d0 L(:, :) = imcompleteCholosky(a) call showArray(L) !do i=2, n ! ! i < kの場合 ! do j=1, i ! if( abs(a(i,j)) < dble(1.0e-10) )then ! cycle ! else ! lld = a(i,j) ! do k=1, j ! lld = lld - L(i,k)*L(j,k)*d(k) ! enddo ! if(d(j)==0.0d0)then ! stop "Error :: d(j)==0.0d0" ! endif ! L(i,j) = 1.0d0/d(j)*lld ! endif ! ld = a(i,i) ! do k=1, i ! ld = ld - L(i,k)*L(i,k)*d(k) ! enddo ! d(i) = ld ! L(i,i) = 1.0d0 ! enddo ! !enddo ! print *, d ! <<<<<<<<<< eps = dble(1.00e-14) er0 = input(default=eps, option=er) r(:) = b - matmul(a, x) call icres(L, d, r, p, n) c1 = dot_product(r, r) init_rr = c1 if (c1 < er0) return beta = 0.0d0 !p(:) = 0.0d0 r0(:) = r(:) w(:) = 0.0d0 u(:) = 0.0d0 t(:) = 0.0d0 t_(:) = 0.0d0 z(:) = 0.0d0 itrmax_ = input(default=100000, option=itrmax) itrmax_ = 10 do itr = 1, itrmax_ call icres(L, d, r, r_k, n) !p(:) = r(:)+beta*(p(:)-u(:) ) ! triple checked. !r0rk=dot_product(r0,r) ap = matmul(a, p) ! triple checked. ap=v, alp = dot_product(r, r_k)/dot_product(p, ap)! triple checked. x(:) = x(:) + alp*ap(:) rk(:) = r(:) r(:) = r(:) - alp*ap(:) !y(:) = t(:) - r(:) - alp*w(:) + alp*ap(:) !triple checked. !t(:) = r(:) - alp*ap(:) ! triple checked. !q(:) = matmul(a,t) ! triple checked. s=q=at=c !if(itr==1)then ! gzi = dot_product(q,t)/dot_product(q,q)! double checked. ! nu = 0.0d0 ! double checked. !else ! val1 = dot_product(y,y)*dot_product(q,t) - dot_product(y,t)*dot_product(q,y) ! val2 = dot_product(q,q)*dot_product(y,y) - dot_product(y,q)*dot_product(q,y) ! ! if( val2==0.0d0 ) stop "Bicgstab devide by zero" ! gzi = val1/val2 ! val1 = dot_product(q,q)*dot_product(y,t) - dot_product(y,q)*dot_product(q,t) ! nu = val1/val2 !triple checked. !endif !u(:) = gzi*ap + nu*(t_(:) -r(:) + beta*u(:) ) !double checked. !z(:) = gzi*r(:) + nu*z(:) - alp*u(:)!double checked. !x(:) = x(:) +alp*p(:) +z(:) !double checked. !r(:) = t(:) -nu*y(:) -gzi*q(:)!double checked. !t_(:)=t(:) rr = dot_product(r, r) print *, rr !print *, 'itr, er =', itr,rr,sqrt(rr)/sqrt(init_rr) !print *, "ans = ",x(:) if (sqrt(rr)/sqrt(init_rr) < er0) exit call icres(L, d, r, r_, n) beta = dot_product(r, r_)/dot_product(rk, r_k) p(:) = r_(:) + beta*p(:) !beta=alp/gzi*dot_product(r0,r)/r0rk !w(:) = q(:) + beta * ap(:) if (itr == itrmax_) then print *, "ERROR :: preconditioned-CG did not converge." end if end do end subroutine !=============================================================== subroutine icres(L, d, r, u, n) real(real64) :: y(n), rly, lu real(real64), intent(inout) :: L(1:n, 1:n), d(1:n), r(1:n), u(1:n) integer(int32), intent(in) :: n integer(int32) :: i, j do i = 1, n rly = r(i) do j = 1, i rly = rly - L(i, j)*y(j) end do y(i) = rly/L(i, i) end do do i = n, 1, -1 lu = 0.0d0 do j = i + 1, n lu = lu + L(j, i)*u(j) end do u(i) = y(i) - d(i)*lu end do end subroutine icres !========================================================== function eigen_3d(tensor) result(eigenvector) real(real64), intent(in) :: tensor(3, 3) integer(int32) :: i, j, n real(real64) :: eigenvector(3, 3), a, b, c, d real(real64) :: eigenvalue(3), mat(3, 3), ev(3), zero(3), unitmat(3, 3) zero(:) = 0.0d0 unitmat(:, :) = 0.0d0 unitmat(1, 1) = 1.0d0 unitmat(2, 2) = 1.0d0 unitmat(3, 3) = 1.0d0 ! get eigen values ! (1) solve cubic equation !https://keisan.casio.jp/exec/user/1305724050 a = 1.0d0 b = -tensor(1, 1) - tensor(2, 2) - tensor(3, 3) c = tensor(1, 1)*tensor(2, 2) + tensor(2, 2)*tensor(3, 3) + tensor(3, 3)*tensor(1, 1) & - tensor(1, 2)*tensor(2, 1) - tensor(2, 3)*tensor(3, 2) - tensor(3, 1)*tensor(1, 3) d = -tensor(1, 1)*tensor(2, 2)*tensor(3, 3) - tensor(1, 2)*tensor(2, 3)*tensor(3, 1) & - tensor(1, 3)*tensor(2, 1)*tensor(3, 2) + tensor(1, 3)*tensor(2, 2)*tensor(3, 1) & + tensor(2, 3)*tensor(3, 2)*tensor(1, 1) + tensor(3, 3)*tensor(1, 2)*tensor(2, 1) eigenvalue = cubic_equation(a, b, c, d) do i = 1, 3 mat(:, :) = eigenvalue(i)*unitmat(:, :) - tensor(:, :) call gauss_jordan_pv(mat, ev, zero, 3) eigenvector(:, 3) = ev(:) end do end function !==================================================================================== subroutine showLinearSolver(obj) class(LinearSolver_), intent(in) :: obj real(real64), allocatable :: A_ij(:, :) integer(int32) :: i, j, n, m n = maxval(obj%Index_I) m = maxval(obj%Index_J) n = maxval((/n, m/)) A_ij = zeros(n, n) do i = 1, size(obj%val) A_ij(obj%Index_I(i), obj%Index_J(i)) = A_ij(obj%Index_I(i), obj%Index_J(i)) + obj%val(i) end do call print(A_ij) end subroutine !==================================================================================== !==================================================================================== function globalMatrixLinearSolver(obj) result(ret) class(LinearSolver_), intent(in) :: obj real(real64), allocatable :: ret(:, :) integer(int32) :: i, j, n, m, p1, p2, domain_id, offset if (allocated(obj%NumberOfNode)) then n = sum(obj%NumberOfNode)*obj%DOF ret = zeros(n, n) do i = 1, size(obj%val) domain_id = obj%row_domain_id(i) offset = sum(obj%NumberOfNode(:domain_id - 1)) !print *, offset p1 = offset*obj%DOF + obj%Index_I(i) domain_id = obj%column_domain_id(i) offset = sum(obj%NumberOfNode(:domain_id - 1)) p2 = offset*obj%DOF + obj%Index_J(i) ret(p1, p2) = ret(p1, p2) + obj%val(i) end do return else n = maxval(obj%Index_I) m = maxval(obj%Index_J) n = maxval((/n, m/)) ret = zeros(n, n) do i = 1, size(obj%val) ret(obj%Index_I(i), obj%Index_J(i)) = ret(obj%Index_I(i), obj%Index_J(i)) + obj%val(i) end do return end if end function !==================================================================================== !==================================================================================== function globalVectorLinearSolver(obj) result(ret) class(LinearSolver_), intent(in) :: obj real(real64), allocatable :: ret(:) integer(int32) :: i, j, n, m, p1, p2, domain_id, offset if (allocated(obj%NumberOfNode)) then n = sum(obj%NumberOfNode)*obj%DOF allocate (ret(n)) ret(:) = 0.0d0 do i = 1, size(obj%b_Index_J) domain_id = obj%b_domain_id(i) offset = sum(obj%NumberOfNode(:domain_id - 1)) p1 = offset*obj%DOF + obj%b_Index_J(i) ret(p1) = ret(p1) + obj%b(i) end do return else ret = obj%b end if end function !==================================================================================== subroutine convertCOOtoCRSLinearSolver(obj, OpenMP) class(LinearSolver_), intent(inout) :: obj integer(int32), allocatable :: buf(:) integer(int32) :: n, nnz, i, nrhs logical, optional, intent(in) :: OpenMP logical :: omp_swich !real(real64),allocatable :: CRS_val(:) !integer(int32),allocatable :: CRS_index_I(:) !integer(int32),allocatable :: CRS_index_J(:) !integer(int32),allocatable :: CRS_row_domain_id(:) !integer(int32),allocatable :: CRS_column_domain_id(:) omp_swich = input(default=.false., option=OpenMP) ! Notice :: COO format data should be created and sorted. ! Further, multi-domain is not supported. nrhs = size(obj%b) + 1 obj%CRS_index_I = int(zeros(nrhs)) if (omp_swich) then !$OMP parallel do private(i) do i = 1, size(obj%index_I) obj%CRS_index_I(obj%index_I(i)) = obj%CRS_index_I(obj%index_I(i)) + 1 end do !$OMP end parallel do buf = obj%CRS_index_I !$OMP parallel do private(i) do i = nrhs - 1, 2, -1 obj%CRS_index_I(i) = sum(buf(1:i - 1)) end do !$OMP end parallel do obj%CRS_index_I(1) = 0 obj%CRS_index_I(nrhs) = size(obj%val) else do i = 1, size(obj%index_I) obj%CRS_index_I(obj%index_I(i)) = obj%CRS_index_I(obj%index_I(i)) + 1 end do buf = obj%CRS_index_I do i = nrhs - 1, 2, -1 obj%CRS_index_I(i) = sum(buf(1:i - 1)) end do obj%CRS_index_I(1) = 0 obj%CRS_index_I(nrhs) = size(obj%val) end if end subroutine ! ####################################################### function matmulCRSLinearSolver(obj, openMP) result(mm) class(LinearSolver_), intent(inout) :: obj real(real64), allocatable :: mm(:) integer(int32) :: i, j logical, optional, intent(in) :: openMP logical :: omp_swich omp_swich = input(default=.false., option=OpenMP) mm = zeros(size(obj%b)) ! Notice :: CRS format data should be created and sorted. if (omp_swich) then !$OMP parallel do private(i) do i = 1, size(obj%b) do j = obj%CRS_Index_I(i) + 1, obj%CRS_Index_I(i + 1) mm(i) = mm(i) + obj%b(obj%Index_J(j))*obj%val(j) end do end do !$OMP end parallel do else do i = 1, size(obj%b) do j = obj%CRS_Index_I(i) + 1, obj%CRS_Index_I(i + 1) mm(i) = mm(i) + obj%b(obj%Index_J(j))*obj%val(j) end do end do end if end function ! ####################################################### function matmulCOOLinearSolver(obj, OpenMP) result(mm) class(LinearSolver_), intent(inout) :: obj real(real64), allocatable :: mm(:) integer(int32) :: i, j logical, optional, intent(in) :: openMP logical :: omp_swich omp_swich = input(default=.false., option=OpenMP) mm = zeros(size(obj%b)) ! Notice :: CRS format data should be created and sorted. if (omp_swich) then !$OMP parallel do private(i) do i = 1, size(obj%val) mm(obj%index_I(i)) = mm(obj%index_I(i)) + obj%val(i)*obj%b(obj%index_J(i)) end do !$OMP end parallel do else do i = 1, size(obj%val) mm(obj%index_I(i)) = mm(obj%index_I(i)) + obj%val(i)*obj%b(obj%index_J(i)) end do end if end function ! ####################################################### ! ####################################################### subroutine getCOOFormatLinearSolver(obj) class(LinearSolver_), intent(inout) :: obj integer(int32) :: i, j, n if (.not. allocated(obj%val)) then if (allocated(obj%a) .and. allocated(obj%b)) then ! it only has obj%a and obj%b ! so, convert it to COO format ! count non-zero values n = 0 do j = 1, size(obj%a, 2) do i = 1, size(obj%a, 1) if (obj%a(i, j) /= 0.0d0) then n = n + 1 end if end do end do obj%val = zeros(n) obj%index_I = int(zeros(n)) obj%index_J = int(zeros(n)) obj%row_domain_id = int(zeros(n)) obj%column_domain_id = int(zeros(n)) obj%b_Index_J = int(zeros(size(obj%b))) obj%b_Domain_ID = int(zeros(size(obj%b))) obj%row_domain_id(:) = 1 obj%column_domain_id(:) = 1 obj%b_domain_id(:) = 1 do i = 1, size(obj%b) obj%b_Index_J(i) = i end do n = 0 do j = 1, size(obj%a, 2) do i = 1, size(obj%a, 1) if (obj%a(i, j) /= 0.0d0) then n = n + 1 obj%val(n) = obj%a(i, j) obj%index_I(n) = i obj%index_J(n) = j end if end do end do end if end if end subroutine ! ####################################################### ! ####################################################### subroutine exportAsCOOLinearSolver(obj, name) class(LinearSolver_), intent(inout) :: obj character(*), intent(in) :: name type(IO_) :: f call f%open(name, "w") call f%write(obj%Index_I, obj%Index_J, obj%val) call f%close() end subroutine ! ####################################################### ! ####################################################### subroutine exportRHSLinearSolver(obj, name) class(LinearSolver_), intent(inout) :: obj character(*), intent(in) :: name type(IO_) :: f call f%open(name, "w") call f%write(obj%b) call f%close() end subroutine ! ####################################################### ! ####################################################### ! FITTING by Stochastic Gradient Descend function fit(f, training_data, params, eta, error, max_itr, use_ratio, logfile, algorithm) result(fit_params) interface function f(x, params) result(ret) use iso_fortran_env real(real64), intent(in) :: x real(real64), intent(in) :: params(:) real(real64) :: ret end function end interface real(real64), intent(in) :: training_data(:, :) !(sample, {input=1, output=2}) real(real64), intent(in) :: params(:), eta real(real64), optional, intent(inout) :: error real(real64), optional, intent(in) :: use_ratio real(real64), allocatable :: grad_params(:), dp(:) integer(int32), optional, intent(in) :: max_itr character(*), optional, intent(in) :: logfile real(real64), allocatable :: fit_params(:) real(real64) :: error_function, error_function_b, error_function_f, err_0, h integer(int32) :: i, j, n, trial, shuffle_id, param_id real(real64) :: eps_val = dble(1.0e-4) real(real64) :: tol = dble(1.0e-9) real(real64) :: sgd_use_ratio character(*), optional, intent(in) :: algorithm character(:), allocatable :: algorithm_type type(IO_) :: log_file real(real64) :: eta_zero real(real64) :: eta_tr integer(int32):: max_trial = 1000000 integer(int32) :: num_select ! type(Random_) :: random integer(int32), allocatable :: selected_data_ids(:) if (present(algorithm)) then algorithm_type = algorithm else algorithm_type = "SGD" end if if (algorithm_type == "SGD") then sgd_use_ratio = 0.050d0 if (present(use_ratio)) then sgd_use_ratio = use_ratio end if num_select = int(dble(size(training_data, 1))*sgd_use_ratio) + 1 selected_data_ids = zeros(num_select) if (present(max_itr)) then max_trial = max_itr end if if (present(logfile)) then call log_file%open(logfile, "w") call log_file%write("# itr error-norm") end if eta_zero = eta fit_params = params error_function = 0.0d0 do i = 1, size(training_data, 1) error_function = error_function + (f(x=training_data(i, 1), params=fit_params) - training_data(i, 2))**2 end do ! compute grad for each params ! by Stochastic Gradient Descend do trial = 1, max_trial error_function = 0.0d0 !$OMP parallel default(shared) !$OMP do reduction(+:error_function) do i = 1, size(training_data, 1) error_function = error_function + (f(x=training_data(i, 1), params=fit_params) - training_data(i, 2))**2 end do !$OMP end do !$OMP end parallel if (trial == 1) then err_0 = error_function end if if (present(error)) then error = error_function/err_0 end if if (error_function/err_0 <= tol) then if (present(logfile)) then call log_file%close() end if exit end if if (trial == max_trial) then if (present(logfile)) then call log_file%close() end if !print *, "ERROR: fit() > SGD did not converge. " return end if ! Robbins-Monro Method eta_tr = eta_zero/dble(trial) ! shuffle training data and select 1 shuffle_id = random%randint(from=1, to=size(training_data, 1)) do i = 1, num_select selected_data_ids(i) = random%randint(from=1, to=size(training_data, 1)) end do ! compute grad grad_params = zeros(size(params)) !$OMP parallel private(shuffle_id,dp,error_function_f,error_function_b) !$OMP do reduction(+:grad_params) do j = 1, num_select shuffle_id = selected_data_ids(j) dp = zeros(size(params)) do param_id = 1, size(params) dp(:) = 0.0d0 dp(param_id) = eps_val error_function_f = 0.0d0 do i = 1, size(training_data, 1) error_function_f = error_function_f + (f(x=training_data(shuffle_id, 1), params=fit_params + dp) & - training_data(shuffle_id, 2))**2 end do dp(param_id) = -eps_val error_function_b = 0.0d0 do i = 1, size(training_data, 1) error_function_b = error_function_b + (f(x=training_data(shuffle_id, 1), params=fit_params + dp) & - training_data(shuffle_id, 2))**2 end do ! numerical derivative grad_params(param_id) = grad_params(param_id) + & (error_function_f - error_function_b)/(2.0d0*eps_val) end do end do !$OMP end do !$OMP end parallel ! average_gradient grad_params = grad_params/dble(num_select) !grad_params = grad_params/norm(grad_params) fit_params = fit_params - eta_tr*grad_params if (norm(grad_params) <= tol) exit grad_params = 0.0d0 if (present(logfile)) then write (log_file%fh, *) trial, error end if end do if (present(logfile)) then call log_file%close() end if else print *, "[ERROR] FIT > no algorithm was found:: ", algorithm_type end if end function function LOBPCG_sparse(A_val, A_col, A_rowptr, lambda_min, tolerance) result(eigen_vectors) real(real64), intent(in) :: A_val(:) integer(int32), intent(in)::A_col(:) integer(int64), intent(in)::A_rowptr(:) real(real64), allocatable :: eigen_vectors(:, :) real(real64), allocatable :: V(:, :), A_(:, :), xtemp(:, :), & lambda_and_x(:, :), lambda(:), bmin(:), x(:, :), r(:, :), x0s(:, :), p(:, :), w(:), & alpha, beta, gamma, XX(:), WW(:), PP(:), SA(:, :), SB(:, :), w_x_p(:, :), WW_XX_PP(:, :), & SB_inv(:, :), x_(:, :), lambda_mat(:, :), lambda_ids(:), Bm(:, :), residual(:), norms(:) integer(int32) :: num_eigen real(real64), optional, intent(in) :: tolerance real(real64) :: tol = dble(1.0e-6) real(real64) :: mu, normval real(real64), intent(inout) :: lambda_min(:) type(Random_) :: random integer(int32) :: i, j, n, id, itr integer(int32) :: m tol = input(default=dble(1.0e-14), option=tolerance) num_eigen = size(lambda_min, 1) n = size(A_rowptr, 1) - 1 if (num_eigen*3 >= n) then print *, "ERROR :: num_eigen*3 should be < n" allocate (eigen_vectors(0, 0)) return end if !m = 3*num_eigen residual = zeros(num_eigen) !https://www.researchgate.net/figure/Algorithm-of-LOBPCG-method-for-matrix-A-Here-the-matrix-T-is-a-preconditioner_fig1_323863889 ! !x0s = random%randn(n,3) !initial guess x0 lambda_min = 0.0d0 x = zeros(n, num_eigen) x = (random%randn(n, num_eigen) - 0.50d0)*2.0d0 x = x*10.0d0 do i = 1, size(x, 1) do j = 1, size(x, 2) if (abs(x(i, j)) < 1.0d0) then x(i, j) = x(i, j)/abs(x(i, j)) end if end do end do itr = 0 do itr = itr + 1 if (itr == 1) then ! step 1 to make initial values if (num_eigen /= 1) then xtemp = X call gram_real(n, num_eigen, Xtemp, X) ! [ok] X(:,:) :: ok! A_ = matmul(transpose(X), crs_matmul(A_val, A_col, A_rowptr, X)) A_ = 0.50d0*A_ + 0.50d0*transpose(A_) ! [ok] A_(:,:) has correct size call eigenValueAndVector(A=A_, lambda=lambda, x=x_, tol=tol*dble(1.0e-3)) lambda_mat = zeros(num_eigen, num_eigen) lambda_mat(:, :) = 0.0d0 do i = 1, size(lambda_mat, 1) lambda_mat(i, i) = lambda(i) end do ! [ok] lambda_mat ! [ok] A_ X = X \Lambda ! X1 = X0 B x = matmul(x, x_) ![ok] size and content of R R = crs_matmul(A_val, A_col, A_rowptr, X) - matmul(X, lambda_mat) else mu = dot_product(x(:, 1), crs_matvec(A_val, A_col, A_rowptr, x(:, 1)))/dot_product(x(:, 1), x(:, 1)) R = crs_matmul(A_val, A_col, A_rowptr, X) - mu*X end if ! 2m次元固有値問題 V = zeros(n, 2*num_eigen) V(:, 1:num_eigen) = X V(:, num_eigen + 1:) = R ! [ok]直交化 xtemp = V call gram_real(n, 2*num_eigen, xtemp, V) X = V(:, 1:num_eigen) R = V(:, num_eigen + 1:) A_ = matmul(transpose(V), crs_matmul(A_val, A_col, A_rowptr, V)) ! try A_ = 0.50d0*A_ + 0.50d0*transpose(A_) ![ok] size(A_,1) & size(A_,2) !2m 次元固有値問題 call eigenValueAndVector(A=A_, lambda=lambda, x=x_, tol=dble(1.0e-14)) ![ok] 下から m 個の固有値と固有ベクトルからなる行列: lambda_ids = linspace([1.0d0, dble(2*num_eigen)], 2*num_eigen) call heapsort(n=2*num_eigen, array=lambda, val=lambda_ids) ![ok] sort lambda_mat = zeros(num_eigen, num_eigen) do i = 1, num_eigen lambda_mat(i, i) = lambda(i) end do ! [ok] check eigen values Bm = zeros(2*num_eigen, num_eigen) do i = 1, num_eigen Bm(:, i) = x_(:, int(lambda_ids(i))) end do ! X2, R2, P2 ! V = {X, R} X = matmul(V, Bm) !V = zeros(n,num_eigen*2 ) V(:, 1:num_eigen) = 0.0d0 V(:, num_eigen + 1:2*num_eigen) = r(:, :) ! R=R1 orR=R2? ! 今はR1と仮定 P = matmul(V, Bm) R = crs_matmul(A_val, A_col, A_rowptr, X) - matmul(X, lambda_mat) !V(:,1:num_eigen) =0.0d0 !V(:,num_eigen+1:2*num_eigen) = r(:,:) ! R=R1 orR=R2? ! 今はR2と仮定 !P = matmul(V,Bm) do i = 1, num_eigen Residual(i) = norm(R(:, i)) end do call print("residual: "+str(maxval(Residual))) cycle else V = zeros(n, 3*num_eigen) V(:, 1:num_eigen) = x(:, :) V(:, num_eigen + 1:2*num_eigen) = r(:, :) V(:, 2*num_eigen + 1:) = p(:, :) xtemp = V call gram_real(n, 3*num_eigen, xtemp, V) x(:, :) = V(:, 1:num_eigen) r(:, :) = V(:, num_eigen + 1:2*num_eigen) p(:, :) = V(:, 2*num_eigen + 1:) A_ = matmul(transpose(V), crs_matmul(A_val, A_col, A_rowptr, V))! try A_ = 0.50d0*A_ + 0.50d0*transpose(A_) x_ = zeros(3*num_eigen, 3*num_eigen) call eigenValueAndVector(A=A_, lambda=lambda, x=x_, tol=dble(1.0e-14)) lambda_ids = linspace([1.0d0, dble(3*num_eigen)], 3*num_eigen) lambda_mat = zeros(num_eigen, num_eigen) call heapsort(n=3*num_eigen, array=lambda, val=lambda_ids) ![ok] sort of lambda do i = 1, num_eigen lambda_mat(i, i) = lambda(i) end do Bm = zeros(3*num_eigen, num_eigen) do i = 1, num_eigen Bm(:, i) = x_(:, int(lambda_ids(i))) end do X = matmul(V, Bm) !V = zeros(n,3*num_eigen ) V(:, 1:num_eigen) = 0.0d0 P = matmul(V, Bm) !V(:,num_eigen+1:2*num_eigen) = r(:,:) !V(:,2*num_eigen+1:) = p(:,:) R = crs_matmul(A_val, A_col, A_rowptr, X) - matmul(X, Lambda_mat) do i = 1, num_eigen Residual(i) = norm(R(:, i)) end do print *, "info: ", itr, maxval(residual), minval(residual) if (maxval(Residual) < tol) then exit end if end if !あとは収束判定 if (itr == 20000) then print *, "ERROR :: LOBPCG >> did not converge!" exit end if end do print *, "residual:", maxval(residual) print *, "itr=", itr do i = 1, num_eigen lambda_min(i) = lambda(i) end do eigen_vectors = x end function !! ################################################################### ! ################################################################### function LOBPCG_dense(A, B, lambda_min) result(eigen_vectors) real(real64), intent(in) :: A(:, :), B(:, :) real(real64), allocatable :: eigen_vectors(:, :) real(real64), allocatable :: V(:, :), A_(:, :), xtemp(:, :), & lambda_and_x(:, :), lambda(:), bmin(:), x(:, :), r(:, :), x0s(:, :), p(:, :), w(:), & alpha, beta, gamma, XX(:), WW(:), PP(:), SA(:, :), SB(:, :), w_x_p(:, :), WW_XX_PP(:, :), & SB_inv(:, :), x_(:, :), lambda_mat(:, :), lambda_ids(:), Bm(:, :), residual(:), norms(:) integer(int32) :: num_eigen real(real64) :: tol = dble(1.0e-14) real(real64) :: mu, normval real(real64), intent(inout) :: lambda_min(:) type(Random_) :: random integer(int32) :: i, j, n, id, itr integer(int32) :: m num_eigen = size(lambda_min, 1) n = size(A, 1) if (num_eigen*3 >= n) then print *, "ERROR :: num_eigen*3 should be < n" allocate (eigen_vectors(0, 0)) return end if !m = 3*num_eigen residual = zeros(num_eigen) !https://www.researchgate.net/figure/Algorithm-of-LOBPCG-method-for-matrix-A-Here-the-matrix-T-is-a-preconditioner_fig1_323863889 ! !x0s = random%randn(n,3) !initial guess x0 lambda_min = 0.0d0 n = size(A, 1) x = zeros(n, num_eigen) x = (random%randn(n, num_eigen) - 0.50d0)*2.0d0 x = x*10.0d0 do i = 1, size(x, 1) do j = 1, size(x, 2) if (abs(x(i, j)) < 1.0d0) then x(i, j) = x(i, j)/abs(x(i, j)) end if end do end do itr = 0 do itr = itr + 1 if (itr == 1) then ! step 1 to make initial values if (num_eigen /= 1) then xtemp = X call gram_real(n, num_eigen, Xtemp, X) ! [ok] X(:,:) :: ok! A_ = matmul(transpose(X), matmul(A, X)) A_ = 0.50d0*A_ + 0.50d0*transpose(A_) ! [ok] A_(:,:) has correct size call eigenValueAndVector(A=A_, lambda=lambda, x=x_, tol=tol) lambda_mat = zeros(num_eigen, num_eigen) lambda_mat(:, :) = 0.0d0 do i = 1, size(lambda_mat, 1) lambda_mat(i, i) = lambda(i) end do ! [ok] lambda_mat ! [ok] A_ X = X \Lambda ! X1 = X0 B x = matmul(x, x_) ![ok] size and content of R R = matmul(A, X) - matmul(X, lambda_mat) else mu = dot_product(x(:, 1), matmul(A(:, :), x(:, 1)))/dot_product(x(:, 1), x(:, 1)) R = matmul(A, X) - mu*X end if ! 2m次元固有値問題 V = zeros(n, 2*num_eigen) V(:, 1:num_eigen) = X V(:, num_eigen + 1:) = R ! [ok]直交化 xtemp = V call gram_real(n, 2*num_eigen, xtemp, V) X = V(:, 1:num_eigen) R = V(:, num_eigen + 1:) A_ = matmul(transpose(V), matmul(A, V)) ! try A_ = 0.50d0*A_ + 0.50d0*transpose(A_) ![ok] size(A_,1) & size(A_,2) !2m 次元固有値問題 call eigenValueAndVector(A=A_, lambda=lambda, x=x_, tol=tol) ![ok] 下から m 個の固有値と固有ベクトルからなる行列: lambda_ids = linspace([1.0d0, dble(2*num_eigen)], 2*num_eigen) call heapsort(n=2*num_eigen, array=lambda, val=lambda_ids) ![ok] sort lambda_mat = zeros(num_eigen, num_eigen) do i = 1, num_eigen lambda_mat(i, i) = lambda(i) end do ! [ok] check eigen values Bm = zeros(2*num_eigen, num_eigen) do i = 1, num_eigen Bm(:, i) = x_(:, int(lambda_ids(i))) end do ! X2, R2, P2 ! V = {X, R} X = matmul(V, Bm) !V = zeros(n,num_eigen*2 ) V(:, 1:num_eigen) = 0.0d0 V(:, num_eigen + 1:2*num_eigen) = r(:, :) ! R=R1 orR=R2? ! 今はR1と仮定 P = matmul(V, Bm) R = matmul(A, X) - matmul(X, lambda_mat) !V(:,1:num_eigen) =0.0d0 !V(:,num_eigen+1:2*num_eigen) = r(:,:) ! R=R1 orR=R2? ! 今はR2と仮定 !P = matmul(V,Bm) do i = 1, num_eigen Residual(i) = norm(R(:, i)) end do call print("residual: "+str(maxval(Residual))) cycle else V = zeros(n, 3*num_eigen) V(:, 1:num_eigen) = x(:, :) V(:, num_eigen + 1:2*num_eigen) = r(:, :) V(:, 2*num_eigen + 1:) = p(:, :) xtemp = V call gram_real(n, 3*num_eigen, xtemp, V) x(:, :) = V(:, 1:num_eigen) r(:, :) = V(:, num_eigen + 1:2*num_eigen) p(:, :) = V(:, 2*num_eigen + 1:) A_ = matmul(transpose(V), matmul(A, V))! try A_ = 0.50d0*A_ + 0.50d0*transpose(A_) x_ = zeros(3*num_eigen, 3*num_eigen) call eigenValueAndVector(A=A_, lambda=lambda, x=x_, tol=tol) lambda_ids = linspace([1.0d0, dble(3*num_eigen)], 3*num_eigen) lambda_mat = zeros(num_eigen, num_eigen) call heapsort(n=3*num_eigen, array=lambda, val=lambda_ids) ![ok] sort of lambda do i = 1, num_eigen lambda_mat(i, i) = lambda(i) end do Bm = zeros(3*num_eigen, num_eigen) do i = 1, num_eigen Bm(:, i) = x_(:, int(lambda_ids(i))) end do X = matmul(V, Bm) !V = zeros(n,3*num_eigen ) V(:, 1:num_eigen) = 0.0d0 P = matmul(V, Bm) !V(:,num_eigen+1:2*num_eigen) = r(:,:) !V(:,2*num_eigen+1:) = p(:,:) R = matmul(A, X) - matmul(X, Lambda_mat) do i = 1, num_eigen Residual(i) = norm(R(:, i)) end do if (maxval(Residual) < tol) then exit end if end if !あとは収束判定 if (itr == 20000) then print *, "ERROR :: LOBPCG >> did not converge!" exit end if end do print *, "residual:", maxval(residual) print *, "itr=", itr do i = 1, num_eigen lambda_min(i) = lambda(i) end do eigen_vectors = x end function !! ################################################################### function LOBPCG_dense_single(A, B, lambda_min) result(eigen_vector) real(real64), intent(in) :: A(:, :), B(:, :) real(real64), allocatable :: eigen_vector(:) real(real64), allocatable :: x0(:), r0(:), p0(:), V(:, :), A_(:, :), xtemp(:, :), & lambda_and_x(:, :), lambda(:), bmin(:), x(:), r(:), x0s(:, :), p(:), w(:), & alpha, beta, gamma, XX(:), WW(:), PP(:), SA(:, :), SB(:, :), w_x_p(:, :), WW_XX_PP(:, :), & SB_inv(:, :), x_(:, :) real(real64) :: tol = dble(1.0e-10) real(real64) :: mu real(real64), intent(inout) :: lambda_min type(Random_) :: random integer(int32) :: i, j, n, id, itr integer(int32) :: m = 3 !https://www.researchgate.net/figure/Algorithm-of-LOBPCG-method-for-matrix-A-Here-the-matrix-T-is-a-preconditioner_fig1_323863889 ! !x0s = random%randn(n,3) !initial guess x0 n = size(A, 1) lambda_min = 0.0d0 n = size(A, 1) xtemp = zeros(n, 3) eigen_vector = zeros(n) !xtemp = x0s x = zeros(n) do i = 1, n x(i) = 1.0d0 end do x = (random%randn(n) - 0.50d0)*2.0d0 x = x/norm(x) mu = dot_product(x, matmul(A, x))/dot_product(x, x) r = matmul(A, x) - mu*x !r = r/norm(r) p = zeros(n) V = zeros(n, 3) V(:, 1) = x(:) V(:, 2) = r(:) V(:, 3) = p(:) xtemp = V call gram_real(n, 2, xtemp, v) x(:) = V(:, 1) r(:) = V(:, 2) V = zeros(n, 3) V(:, 1) = x(:) V(:, 2) = r(:) V(:, 3) = p(:) ! itr = 0 do itr = itr + 1 A_ = matmul(transpose(V), matmul(A, V)) A_ = 0.50d0*A_ + 0.50d0*transpose(A_) !lambda = eigenValue(A_,tol=tol) call eigenValueAndVector(A=A_, lambda=lambda, x=x_, tol=tol) ! 固有ベクトルを正規化 do i = 1, size(x_, 2) x_(:, i) = x_(:, i)/norm(x_(:, i)) end do if (itr == 1) then id = minvalID(lambda(1:2)) else id = minvalID(lambda(1:3)) end if print *, "---check---" ! 最小固有ベクトル bmin = x_(:, id) x = matmul(V, bmin) mu = lambda(id) r = matmul(A, x) - mu*x(:) p = zeros(n) p(:) = V(:, 2)*bmin(2) + V(:, 3)*bmin(3) if (norm(r) < tol .and. norm(x) /= 0.0d0) exit print *, "error", norm(r) ! 直交化 V = zeros(n, 3) V(:, 1) = x(:) V(:, 2) = r(:) V(:, 3) = p(:) call gram_real(n, 3, V, xtemp) x(:) = xtemp(:, 1) r(:) = xtemp(:, 2) p(:) = xtemp(:, 3) V(:, 1) = x(:) V(:, 2) = r(:) V(:, 3) = p(:) !print *, dot_product(x,r) !print *, dot_product(r,p) !print *, dot_product(x,p) !x = x/norm(x) !r = r/norm(r) !p = p/norm(p) !call print(xtemp) !stop if (itr == 10000) then print *, "ERROR :: LOBPCG >> did not converge!" exit end if end do print *, "residual:", norm(r) print *, "itr=", itr lambda_min = lambda(id) eigen_vector = x end function !! ################################################################### subroutine gram_real(m, n, mat_v, mat_v_out) ! cited from !http://park.itc.u-tokyo.ac.jp/kato-yusuke-lab/nagai/note_141009_LOBPCG.pdf implicit none ! m : ベクトルの次元 ! n : ベクトルの本数 integer, intent(in)::m, n real(8), intent(in)::mat_v(1:m, 1:n) real(8), intent(out)::mat_v_out(1:m, 1:n) integer::i, j real(8)::v(1:m), nai, vi(1:m), viold(1:m) real(8)::norm mat_v_out = mat_v do i = 1, n viold = mat_v_out(:, i) do j = 1, i - 1 nai = dot_product(mat_v_out(1:m, j), viold(1:m)) vi = viold - nai*mat_v_out(:, j) viold = vi end do norm = sqrt(dble(dot_product(viold, viold))) if (norm == 0.0d0) then ! DEBUG Right? print *, "ERROR gram_real :: norm-zero" mat_v_out(j, j) = 1.0d0 ! !stop !mat_v_out(:,j) = viold else mat_v_out(:, j) = viold/norm end if end do return end subroutine gram_real !function LOBPCG(A,B,num_eigens,err) result(eigens) ! type(FEMSolver_),intent(in) :: A, B ! complex(real64),allocatable :: eigens(:,:) ! integer(int32),intent(in) :: num_eigens ! real(real64),optional,intent(in) :: err ! real(real64) :: error_tolerance ! integer(int32) :: i,j,m,nx,nv ! ! ! Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) ! ! https://arxiv.org/pdf/1704.07458.pdf ! ! ! if(.not.allocated(A%CRS_val) )then ! print *, "[ERROR] eigFEMSolver >> (.not.allocated(this%CRS_val) )" ! return ! endif ! ! ! error_tolerance = input(default=dble(1.0e-14),option=err) ! ! m = size(A%CRS_Index_Row)-1 ! nx= num_eigens ! allocate(eigens(0:num_eigens,1:m)) ! !eigens(0,1:m) :: eigen values ! !eigens(1:nx,1:m) :: eigen vectors from 1st to nx-th ! eigens(:,:) = 0.0d0 ! ! ! ! ! !end function ! ################################################################### !subroutine RayleighRitz(S,C,Theta) ! complex(real64),intent(in) :: S(:,:) ! complex(real64),allocatable,intent(out) :: C(:,:),Theta(:,:) ! complex(real64),allocatable :: D(:,:) ! ! !Rayleigh-Ritz procedure ! D = matmul( transpose(S), ) ! !end subroutine subroutine to_Dense(CRS_val, CRS_col, CRS_rowptr, DenseMatrix) real(real64), allocatable, intent(inout) :: DenseMatrix(:, :) real(real64), allocatable, intent(in) :: CRS_val(:) integer(int32), allocatable, intent(in) :: CRS_col(:) integer(int64), allocatable, intent(in) :: CRS_rowptr(:) integer(int64) :: nonzero_count, i, j, k, n n = size(CRS_rowptr) - 1 allocate (DenseMatrix(n, n)) DenseMatrix(:, :) = 0.0d0 do i = 1, size(CRS_rowptr) - 1 do j = CRS_rowptr(i), CRS_rowptr(i + 1) - 1 DenseMatrix(i, CRS_col(j)) = CRS_val(j) end do end do end subroutine !subroutine to_CRS(DenseMatrix,CRS_val,CRS_col,CRS_rowptr) ! real(real64),intent(in) :: DenseMatrix(:,:) ! real(real64),allocatable,intent(inout) :: CRS_val(:) ! integer(int32),allocatable,intent(inout) :: CRS_col(:) ! integer(int64),allocatable,intent(inout) :: CRS_rowptr(:) ! ! integer(int64) :: nonzero_count,i,j,k,n ! ! nonzero_count = 0 ! do i=1,size(DenseMatrix,1) ! do j=1,size(DenseMatrix,2) ! if(DenseMatrix(i,j)/=0.0d0 ) then ! nonzero_count = nonzero_count + 1 ! endif ! enddo ! enddo ! CRS_val = zeros(nonzero_count) ! CRS_col = int(linspace([1.0d0,dble(nonzero_count)], nonzero_count)) ! CRS_rowptr = int(zeros(size(DenseMatrix,1)+1) ) ! ! nonzero_count = 0 ! do i=1,size(DenseMatrix,1) ! do j=1,size(DenseMatrix,2) ! if(DenseMatrix(i,j)/=0.0d0 ) then ! nonzero_count = nonzero_count + 1 ! CRS_val(nonzero_count) = DenseMatrix(i,j) ! CRS_col(nonzero_count) = j ! CRS_rowptr(i) = CRS_rowptr(i) + 1 ! endif ! enddo ! enddo ! ! n = CRS_rowptr(1) ! !CRS_rowptr(1) =0 ! do i=1,size(CRS_rowptr)-1 ! CRS_rowptr(i+1) = CRS_rowptr(i+1) + CRS_rowptr(i) ! enddo ! do i=size(CRS_rowptr)-1,1,-1 ! CRS_rowptr(i+1) = CRS_rowptr(i) ! enddo ! do i=size(CRS_rowptr),1,-1 ! CRS_rowptr(i) = CRS_rowptr(i) +1 ! enddo ! CRS_rowptr(1) = 1 ! ! !end subroutine ! function crs_matmul_generic(CRS_value, CRS_col, CRS_row_ptr, old_vectors) result(new_vectors) real(real64), intent(in) :: CRS_value(:), Old_vectors(:, :) integeR(int32), intent(in):: CRS_col(:) integeR(int64), intent(in):: CRS_row_ptr(:) real(real64), allocatable :: new_vectors(:, :) integer(int64) :: i, j, n, gid, lid, row, CRS_id, col, m !> x_i = A_ij b_j n = size(CRS_row_ptr) - 1 m = size(old_vectors, 2) if (size(old_vectors, 1) /= n) then print *, "ERROR crs_matmul :: inconsistent size for old_vectors" return end if new_vectors = zeros(n, m) !!$OMP parallel do default(shared) private(CRS_id,col) do row = 1, n do CRS_id = CRS_row_ptr(row), CRS_row_ptr(row + 1) - 1 col = CRS_col(CRS_id) do j = 1, m !!$OMP atomic new_vectors(row, j) = new_vectors(row, j) & + CRS_value(CRS_id)*old_vectors(col, j) end do end do end do !!$OMP end parallel do end function ! ################################################################### function crs_matmul_for_CRStype(CRS, old_vectors) result(new_vectors) type(CRS_), intent(in) :: CRS real(real64), intent(in) :: Old_vectors(:, :) real(real64), allocatable :: new_vectors(:, :) new_vectors = crs_matmul_generic( & CRS_value=CRS%val, & CRS_col=CRS%col_idx, & CRS_row_ptr=CRS%row_ptr, & old_vectors=old_vectors) end function function crs_matvec_generic(CRS_value, CRS_col, CRS_row_ptr, old_vector) result(new_vector) real(real64), intent(in) :: CRS_value(:), Old_vector(:) integeR(int32), intent(in):: CRS_col(:) integeR(int64), intent(in):: CRS_row_ptr(:) real(real64), allocatable :: new_vector(:) integer(int64) :: i, j, n, gid, lid, row, CRS_id, col !> x_i = A_ij b_j n = size(CRS_row_ptr) - 1 if (size(old_vector) /= n) then print *, "ERROR crs_matvec :: inconsistent size for old_vector" return end if new_vector = zeros(n) !$OMP parallel do default(shared) private(CRS_id,col) do row = 1, n do CRS_id = CRS_row_ptr(row), CRS_row_ptr(row + 1) - 1 col = CRS_col(CRS_id) !$OMP atomic new_vector(row) = new_vector(row) + CRS_value(CRS_id)*old_vector(col) end do end do !$OMP end parallel do end function ! ################################################################### subroutine sub_crs_matvec_generic(CRS_value, CRS_col, CRS_row_ptr, old_vector, new_vector, precondition) real(real64), intent(in) :: CRS_value(:), Old_vector(:) integeR(int32), intent(in):: CRS_col(:) integeR(int64), intent(in):: CRS_row_ptr(:) type(CRS_), optional, intent(in) :: precondition real(real64), allocatable, intent(inout) :: new_vector(:) real(real64), allocatable :: precon_old_vector(:) integer(int64) :: i, j, n, gid, lid, row, CRS_id, col !> x_i = A_ij b_j n = size(CRS_row_ptr) - 1 if (size(old_vector) /= n) then print *, "ERROR crs_matvec :: inconsistent size for old_vector" return end if if (.not. allocated(new_vector)) then new_vector = zeros(n) else new_vector(:) = 0.0d0 end if if (present(precondition)) then call precondition%ILU_matvec(old_vector=old_vector, new_vector=precon_old_vector) !$OMP parallel do default(shared) private(CRS_id,col) do row = 1, n do CRS_id = CRS_row_ptr(row), CRS_row_ptr(row + 1) - 1 col = CRS_col(CRS_id) !$OMP atomic new_vector(row) = new_vector(row) + CRS_value(CRS_id)*precon_old_vector(col) end do end do !$OMP end parallel do else !$OMP parallel do default(shared) private(CRS_id,col) do row = 1, n do CRS_id = CRS_row_ptr(row), CRS_row_ptr(row + 1) - 1 col = CRS_col(CRS_id) !$OMP atomic new_vector(row) = new_vector(row) + CRS_value(CRS_id)*old_vector(col) end do end do !$OMP end parallel do end if end subroutine ! ################################################################### function crs_matvec_for_CRStype(CRS, old_vector) result(new_vector) type(CRS_), intent(in) :: CRS real(real64), intent(in) :: Old_vector(:) real(real64), allocatable :: new_vector(:) new_vector = crs_matvec_generic( & CRS_value=CRS%val, & CRS_col=CRS%col_idx, & CRS_row_ptr=CRS%row_ptr, & old_vector=old_vector) end function subroutine sub_crs_matvec_for_CRStype(CRS, old_vector, new_vector, precondition) type(CRS_), intent(in) :: CRS type(CRS_), optional, intent(in) :: precondition real(real64), intent(in) :: Old_vector(:) real(real64), allocatable, intent(inout) :: new_vector(:) call sub_crs_matvec_generic( & CRS_value=CRS%val, & CRS_col=CRS%col_idx, & CRS_row_ptr=CRS%row_ptr, & old_vector=old_vector, & new_vector=new_vector, & precondition=precondition) end subroutine pure function eyesMatrix(rank1, rank2) result(ret) integer(int32), intent(in) ::rank1, rank2 real(real64), allocatable :: ret(:, :) integer(int32) :: i, min_rank allocate (ret(rank1, rank2)) ret(:, :) = 0.0d0 min_rank = minval([rank1, rank2]) do i = 1, min_rank if (rank2 > i) exit ret(i, i) = 1.0d0 end do end function ! ################################################################### subroutine reduce_crs_matrix(CRS_val, CRS_col, CRS_rowptr, remove_IDs) real(real64), allocatable, intent(inout) :: CRS_val(:) integer(int32), allocatable, intent(inout) :: CRS_col(:) integer(int64), allocatable, intent(inout) :: CRS_rowptr(:) real(real64), allocatable :: new_CRS_val(:) integer(int32), allocatable :: new_CRS_col(:) integer(int64), allocatable :: new_CRS_rowptr(:) integer(int32), intent(in) :: remove_IDs(:) integer(int32), allocatable :: remove_IDs_sorted(:) integer(int32) :: i, old_row_id, old_col_id, old_col_ptr, count_m1 integer(int32) :: j integer(int32) :: new_row_id, new_col_id, new_col_ptr, remove_id_ptr, rem_count, k, l, m integer(int32), allocatable :: new_id_from_old_id(:) ! remove_ids should be unique and sorted remove_IDs_sorted = remove_IDs remove_IDs_sorted = unique(remove_IDs_sorted) call heapsort(n=size(remove_IDs_sorted), array=remove_IDs_sorted) ! allocate new_CRS_val = CRS_val new_CRS_col = CRS_col new_CRS_rowptr = CRS_rowptr !new_CRS_val = zeros(size(CRS_val) - size(remove_IDs_sorted) ) !new_CRS_col = zeros(size(CRS_col) - size(remove_IDs_sorted) ) !new_CRS_rowptr = zeros(size(CRS_rowptr) - size(remove_IDs_sorted) ) !! new_id_from_old_id = zeros(size(CRS_rowptr) - 1) do i = 1, size(CRS_rowptr) - 1 new_id_from_old_id(i) = i end do do i = 1, size(remove_ids_sorted) new_id_from_old_id(remove_ids_sorted(i) + 1:) = new_id_from_old_id(remove_ids_sorted(i) + 1:) - 1 new_id_from_old_id(remove_ids_sorted(i)) = -1 end do ! if a id is listed in remove_ids_sorted, ignore ! new_row_id = 0 new_col_id = 0 new_col_ptr = 0 remove_id_ptr = 1 ! create new_CRS_col ! only for column ! blanks are indicated by -1 ! 当該は-1して，前送り do i = 1, size(CRS_rowptr) - 1 do j = CRS_rowptr(i), CRS_rowptr(i + 1) - 1 old_col_ptr = j old_col_id = CRS_col(old_col_ptr) new_CRS_col(old_col_ptr) = new_id_from_old_id(old_col_id) end do end do ! 列に-1 do i = 1, size(remove_IDs_sorted) do j = CRS_rowptr(remove_IDs_sorted(i)), CRS_rowptr(remove_IDs_sorted(i) + 1) - 1 old_col_ptr = j new_CRS_col(old_col_ptr) = -1 end do end do ! renew row_ptr new_CRS_rowptr = 0.0d0 do i = 1, size(CRS_rowptr) - 1 do j = CRS_rowptr(i), CRS_rowptr(i + 1) - 1 if (new_CRS_col(j) /= -1) then new_CRS_rowptr(i) = new_CRS_rowptr(i) + 1 end if end do end do do i = 1, size(new_CRS_rowptr) - 1 new_CRS_rowptr(i + 1) = new_CRS_rowptr(i + 1) + new_CRS_rowptr(i) end do do i = size(new_CRS_rowptr), 2, -1 new_CRS_rowptr(i) = new_CRS_rowptr(i - 1) + 1 end do new_CRS_rowptr(1) = 1 ! set -1 to removable rows do i = 1, size(remove_ids_sorted) new_CRS_rowptr(remove_ids_sorted(i)) = -1 end do ! renew CRS_val k = 0 do i = 1, sizE(new_CRS_col) if (new_CRS_col(i) == -1) then k = k + 1 end if end do new_CRS_val = zeros(size(CRS_val) - k) k = 0 do i = 1, size(new_CRS_col) if (new_CRS_col(i) == -1) then cycle else k = k + 1 new_CRS_val(k) = CRS_val(i) end if end do call searchAndRemove(new_CRS_col, eq=-1) call searchAndRemove(new_CRS_rowptr, eq=-1) deallocate (CRS_val) deallocate (CRS_col) deallocate (CRS_rowptr) CRS_val = new_CRS_val CRS_col = new_CRS_col CRS_rowptr = new_CRS_rowptr end subroutine ! ################################################################### function UpperTriangularMatrix(CRS_val, CRS_col, CRS_rowptr) result(UP) real(real64), intent(in) :: CRS_val(:) integer(int32), intent(in) :: CRS_col(:) integer(int64), intent(in) :: CRS_rowptr(:) integer(int64) :: i, j, col, row, N, offset real(real64) :: val real(real64), allocatable :: UP(:) N = size(CRS_rowptr) - 1 UP = zeros(N*(N + 1)/2) !$OMP parallel default(shared) private(col,val,i) !$OMP do do row = 1, N do i = CRS_rowptr(row), CRS_rowptr(row + 1) - 1 col = CRS_col(i) val = CRS_val(i) if (row <= col) then !!$OMP atomic UP(row + (col - 1)*col/2) = val end if end do end do !$OMP end do !$OMP end parallel end function ! ################################################################### function fillby(element, vec_size, num_repeat) result(new_vec) real(real64), intent(in) :: element(:) integer(int32), optional, intent(in) :: vec_size, num_repeat real(real64), allocatable :: new_vec(:) integer(int32) :: i, j if (present(vec_size)) then allocate (new_vec(vec_size)) i = 0 do do j = 1, 3 i = i + 1 if (i > size(new_vec)) then return end if new_vec(i) = element(j) end do end do elseif (present(num_repeat)) then allocate (new_vec(num_repeat*size(element))) do i = 1, num_repeat new_vec((i - 1)*size(element) + 1:i*size(element)) = element(1:size(element)) end do else new_vec = element end if end function ! ################################################################ function unique(old_vec) result(new_vec) integer(int32), intent(in) :: old_vec(:) integer(int32), allocatable :: new_vec(:), remove_is_one(:) integer(int32) :: n_size, i, id, j n_size = size(old_vec) ! caution; O(N^2) remove_is_one = zeros(n_size) do i = 1, n_size if (remove_is_one(i) == 1) cycle do j = i + 1, n_size if (old_vec(i) == old_vec(j)) then remove_is_one(j) = 1 end if end do end do new_vec = zeros(n_size - sum(remove_is_one)) id = 0 do i = 1, n_size if (remove_is_one(i) == 1) then cycle else id = id + 1 new_vec(id) = old_vec(i) end if end do end function subroutine Lanczos(A, V, T) real(real64), intent(in) :: A(:, :) real(real64), allocatable, intent(inout) :: V(:, :), T(:, :) real(real64), allocatable :: v1(:), w1(:), vj(:), wj(:) real(real64) :: alp1, beta_j, alp_j type(Random_) :: random integer(int32) :: i, N, j, m ! https://en.wikipedia.org/wiki/Lanczos_algorithm N = size(A, 1) V = zeros(N, N) T = zeros(N, N) v1 = 2.0d0*random%randn(N) - 1.0d0 v1 = v1/norm(v1) w1 = matmul(A, v1) alp1 = dot_product(w1, v1) w1 = w1 - alp1*v1 beta_j = norm(w1) T(1, 1) = alp1 V(:, 1) = v1 do i = 2, N beta_j = norm(w1) if (beta_j /= 0.0d0) then vj = w1/beta_j else vj = 2.0d0*random%randn(N) - 1.0d0 stop end if wj = matmul(A, vj) alp_j = dot_product(wj, vj) wj = wj - alp_j*vj - beta_j*v1 v1 = vj w1 = wj T(i, i) = alp_j T(i, i - 1) = beta_j T(i - 1, i) = beta_j V(:, i) = vj end do end subroutine ! ##################################################### subroutine PBiCGSTAB_CRS(CRS, x, b, itrmax, er, relative_er, debug, ILU_MATRIX) type(CRS_), intent(inout) :: CRS integer(int32), intent(inout) :: itrmax real(real64), intent(inout) :: b(:), er real(real64), optional, intent(in) :: relative_er real(real64), intent(inout) :: x(:) logical, optional, intent(in) :: debug type(CRS_), optional, intent(inout) :: ILU_MATRIX call bicgstab_CRS_ILU(a=CRS%val, ptr_i=CRS%row_ptr, index_j=CRS%col_idx, & x=x, b=b, itrmax=itrmax, er=er, relative_er=relative_er, debug=debug, ILU_MATRIX=ILU_MATRIX) end subroutine ! ##################################################### subroutine bicgstab_CRS_ILU(a, ptr_i, index_j, x, b, itrmax, er, relative_er, debug, ILU_MATRIX) integer(int32), intent(inout) :: index_j(:), itrmax integer(int64), intent(inout) :: ptr_i(:) real(real64), intent(inout) :: a(:), b(:), er real(real64), optional, intent(in) :: relative_er real(real64), intent(inout) :: x(:) logical, optional, intent(in) :: debug type(CRS_), optional, intent(inout) :: ILU_MATRIX type(CRS_) :: Pmat logical :: speak = .false. integer(int32) itr, i, j, n real(real64) alp, bet, c1, c2, c3, ev, vv, rr, er0, init_rr, re_er0 real(real64), allocatable:: r(:), r0(:), p(:), y(:), e(:), v(:), pa(:), ax(:), pd(:), ps(:) er0 = er if (present(debug)) then speak = debug end if call Pmat%init(val=a, row_ptr=int(ptr_i, int64), col_idx=index_j) if (speak) print *, "BiCGSTAB ILU(0) STARTED >> " if (present(ILU_MATRIX)) then if (allocated(ILU_MATRIX%val)) then Pmat = ILU_MATRIX else call Pmat%ILU(0, debug=debug) ILU_MATRIX = Pmat end if else call Pmat%ILU(0, debug=debug) end if n = size(b) if (speak) print *, "BiCGSTAB STARTED >> DOF:", n allocate (r(n), r0(n), p(n), y(n), e(n), v(n)) r(:) = b(:) if (speak) print *, "BiCGSTAB >> [1] initialize" call sub_crs_matvec(CRS_value=a, CRS_col=index_j, & CRS_row_ptr=ptr_i, old_vector=x, new_vector=ax) r = b - ax if (speak) print *, "BiCGSTAB >> [2] dp1" c1 = dot_product(r, r) !call omp_dot_product(r,r,c1) init_rr = c1 !if(speak) print *, "BiCGSTAB >> |r|^2 = ",init_rr if (c1 < er0) return p(:) = r(:) r0(:) = r(:) do itr = 1, itrmax if (speak) print *, "BiCGSTAB >> ["//str(itr)//"] initialize" c1 = dot_product(r0, r) !call omp_dot_product(r0,r,c1) y(:) = 0.0d0 ! compute y = A P^{-1} d (p:=d) call sub_crs_matvec(CRS_value=a, CRS_col=index_j, & CRS_row_ptr=ptr_i, old_vector=p, new_vector=y, precondition=Pmat) c2 = dot_product(r0, y) !call omp_dot_product(r0,y,c2) alp = c1/c2 ! compute e := s = r - alpa * A P^{-1} d (p:=d) e(:) = r(:) - alp*y(:) v(:) = 0.0d0 ! v = A P^{-1} e = A P^{-1} s call sub_crs_matvec(CRS_value=a, CRS_col=index_j, & CRS_row_ptr=ptr_i, old_vector=e, new_vector=v, precondition=Pmat) if (speak) print *, "BiCGSTAB >> ["//str(itr)//"] half" ev = dot_product(e, v) vv = dot_product(v, v) !call omp_dot_product(e,v,ev) !call omp_dot_product(v,v,vv) if (vv == 0.0d0) stop "Bicgstab devide by zero" ! w = c3 = c3 = ev/vv !https://cattech-lab.com/science-tools/simulation-lecture-mini-cg/ call Pmat%ILU_matvec(old_vector=p, new_vector=pd) call Pmat%ILU_matvec(old_vector=e, new_vector=ps) x(:) = x(:) + alp*pd(:) + c3*ps(:) !x(:) = x(:) + alp * p(:) + c3 * e(:) r(:) = e(:) - c3*v(:) rr = dot_product(r, r) !call omp_dot_product(r,r,rr) if (itr == 1) then re_er0 = rr end if if (speak) then print *, sqrt(rr) end if if (present(relative_er)) then if (sqrt(rr/re_er0) < relative_er) then exit end if end if ! write(*,*) 'itr, er =', itr,rr if (sqrt(rr) < er0) exit c1 = dot_product(r0, r) !call omp_dot_product(r0,r,c1) bet = c1/(c2*c3) if ((c2*c3) == 0.0d0) stop "Bicgstab devide by zero" p(:) = r(:) + bet*(p(:) - c3*y(:)) end do end subroutine !=============================================================== ! ################################################### function SGD_objective_function(fx, training_data_x, training_data_fx, & params) result(ret) interface function fx(x, params) result(ret) use iso_fortran_env real(real64), intent(in) :: x real(real64), intent(in) :: params(:) real(real64) :: ret end function end interface real(real64), intent(in) :: training_data_x(:) real(real64), intent(in) :: training_data_fx(:) real(real64), intent(in) :: params(:) real(real64) :: ret integer(int32) :: i ret = 0.0d0 do i = 1, size(training_data_x) ret = ret + 1.0d0/dble(size(training_data_x))*(training_data_fx(i) - fx(x=training_data_x(i), params=params))**2 end do end function ! ################################################### function StochasticGradientDescent_scalar_function(fx, & training_data_x, training_data_fx, init_params, eta, max_itr, tol, logfile_idx) result(params) interface function fx(x, params) result(ret) use iso_fortran_env real(real64), intent(in) :: x real(real64), intent(in) :: params(:) real(real64) :: ret end function end interface real(real64), intent(in) :: training_data_x(:) real(real64), intent(in) :: training_data_fx(:) real(real64), intent(in) :: init_params(:) real(real64), intent(in) :: eta, tol integer(int32), intent(in) :: max_itr integer(int32), optional, intent(in) :: logfile_idx real(real64), allocatable :: params(:), d_params(:) integer(int32) :: itr, idx real(real64) :: num integer(int32) :: SGD_log params = init_params if (present(logfile_idx)) then open (newunit=SGD_log, file="SGD"+zfill(logfile_idx, 6) + ".log", status="replace") else open (newunit=SGD_log, file="SGD.log", status="replace") end if do itr = 1, max_itr call random_number(num) idx = int(num*dble(size(training_data_x)) + 1.0d0) !print *, idx d_params = -eta*SGD_gradient( & fx=fx, & training_data_x=training_data_x(idx), & training_data_fx=training_data_fx(idx), & params=params & ) write (SGD_log, *) itr, SGD_objective_function( & fx=fx, & training_data_x=training_data_x, & training_data_fx=training_data_fx, & params=params & ) ! 勾配がNaNであれば無視 if (maxval(d_params(:) - d_params(:)) /= 0.00d0 .or. minval(d_params(:) - d_params(:)) /= 0.00d0) cycle if (maxval(abs(d_params)) < tol) exit params = params + d_params end do close (SGD_log) end function ! ################################################### ! ################################################### function SGD_gradient(fx, training_data_x, training_data_fx, & params) result(grad) interface function fx(x, params) result(ret) use iso_fortran_env real(real64), intent(in) :: x real(real64), intent(in) :: params(:) real(real64) :: ret end function end interface real(real64), intent(in) :: training_data_x real(real64), intent(in) :: training_data_fx real(real64), intent(in) :: params(:) real(real64), allocatable :: grad(:) real(real64), allocatable :: d_params(:) real(real64) :: epsilon integer(int32) :: i, j epsilon = dble(1.0e-4) grad = 0.0d0*params d_params = 0.0d0*params do j = 1, size(grad) d_params = 0.0d0*params d_params(j) = epsilon grad(j) = -2.0d0*(training_data_fx - fx(x=training_data_x, params=params)) & *(fx(x=training_data_x, params=params + d_params) - fx(x=training_data_x, params=params))/epsilon end do end function ! ################################################### subroutine importCRS_LinearSolver(this, CRS) class(LinearSolver_), intent(inout) :: this type(CRS_), intent(in) :: CRS this%CRS_val = CRS%val this%CRS_Index_I = CRS%col_idx this%CRS_Index_J = CRS%row_ptr end subroutine end module