-3
Я написал код Matlab для метода CGR, который отлично работает. Вот оно:Как правильно преобразовать код MATLAB в Fortran 95 или выше?
MATLAB КОД
% Generalized Conjugate Residual Algorithm
% Solves M*x=b (J*dx=-F)
function [x,converged]=gcr_Mfree(F,Cnew,Cold,Fold,alpha);
% define global variables
global h n M;
%Change above
b = -F ;
tol = 5.E-3 ;
maxiter = length(b) ;
x = zeros(maxiter,1) ;
r = b ; % we need this when calculating r_new
normr(1) = norm(r) ;
alpha2 = 1e-6 ; % alpha2 is the epsilon in the report!
for iter=1:maxiter
%Get preliminary search direction
p(:,iter) = r ;
%% Approximate the Jacobian(M) residual product
%alpha2 = 1.E-6 ;
Cnew(1:n) = Cnew(1:n)+alpha2*r ;
% Cnew(1) = Cold(1) ;
F2 = Funct2(Cnew,Cold,Fold,alpha) ;
Mr = 1/alpha2*(F2-F) ;
Mp(:,iter) = Mr ;
%% Orthogonalize search direction
for j=1:iter-1
p(:,iter) = p(:,iter) -Mp(:,j)'*Mp(:,iter)*p(:,j) ;
Mp(:,iter) = Mp(:,iter) -Mp(:,j)'*Mp(:,iter)*Mp(:,j) ;
end
%norm(Mp(:,iter))
%Normalize search direction
p(:,iter) = p(:,iter)/(norm(Mp(:,iter))) ;
Mp(:,iter) = Mp(:,iter)/(norm(Mp(:,iter))) ;
%Update solution and residual
alpha2 = r'*(Mp(:,iter))/((Mp(:,iter))'*(Mp(:,iter))) ;
%alpha2
x = x+alpha2*p(:,iter) ;
r = r-alpha2*(Mp(:,iter)) ;
%Check convergence
normr=norm(r);
%fprintf('norm(r) = %g iter = %g\n',normr,iter+1);
if(normr < tol)
%fprintf('GCR Converged, iter = %g\n',iter+1);
converged=1;
break;
end
end %for loop
if(normr > tol)
fprintf('SOLUTION DID NOT CONVERGE!!');
converged=0;
end
А вот моя последняя версия Fortran кода:
SUBROUTINE gcr_Mfree(F2 ,Cnew,Cold,C_Fold,IGG,JGG,x ,converged)
!**** FUNCTIONS TO BE SOLVED ****
! Generalized Conjugate Residual Algorithm
! Solves M*x=b (J*dx=-F)
! Use somemodule
! Arguments declarations
IMPLICIT REAL*8 (A-H,O-Z)
INTEGER :: IGG,JGG
real*8, dimension(:), ALLOCATABLE :: x
real*8, intent(out) :: converged !<
REAL*8, DIMENSION(:,:), ALLOCATABLE :: F2,F22,Cnew,Cold,C_Fold,p,Mp
! Variable declarations
real*8 :: tol
real*8, DIMENSION(:), ALLOCATABLE :: r,b,Mr,F2V,F22V,CnewV,ColdV,C_FoldV,alpha2
integer :: i,j,maxiter,normr,iter
MASK = SIZE(F2)
ALLOCATE(F2V(1:MASK))
ALLOCATE(CnewV(1:MASK))
ALLOCATE(ColdV(1:MASK))
ALLOCATE(C_FoldV(1:MASK))
!*********** RESHAPING MATRICES TO VECTORS.........
F2V = RESHAPE(F2 ,(/MASK/))
CnewV = RESHAPE(Cnew,(/MASK/))
ColdV = RESHAPE(Cold,(/MASK/))
C_FoldV = RESHAPE(C_Fold,(/MASK/))
b = -F2V !(why minus?)
tol = 5.E-2
alpha2 = 1e-6 ! alpha2 is the epsilon in the report!
maxiter = size(b)
DEALLOCATE(x)
ALLOCATE(x(1:maxiter))
x = 0.0D0 !GUIdE: ""m2f: x = zeros(maxiter,1)""
r = b ! we need this when calculating r_new
! normr(1) = norm2(r) !!!! Norm
DO iter=1,maxiter
!Get preliminary search direction
DO i=1,iter
p(:,i) = r
ENDDO
! Approximate the Jacobian(M) residual product
CnewV = CnewV + alpha2*r
Cnew = RESHAPE(CnewV,(/IGG,JGG/))
F22 = CNF2(Cnew,Cold,Fold,DT)
F22V= RESHAPE(F22V,(/MASK/))
Mr = (1/alpha2)*(F22V-F2V) !GUIDE: (The apporximated Jacobian matrix)
DO i=1,iter
Mp(:,i) = Mr
ENDDO
!! Orthogonalize search direction
do j=1,iter-1
!m2f: p(:,iter) = p(:,iter) - Mp(:,j)'* Mp(:,iter) * p(:,j)
p(:,iter) = p(:,iter) - Mp(j,:) * Mp(:,iter) * p(:,j)
!m2f: Mp(:,iter) = Mp(:,iter) - Mp(:,j)'* Mp(:,iter) * Mp(:,j)
Mp(:, iter) = Mp(:, iter) - sum(Mp(:, j) * Mp(:, iter)) * Mp(:, j)
end do
!norm(Mp(:,iter))
!Normalize search direction
ALLOCATE(p(1:IG,1:JG))
ALLOCATE(Mp(1:IG,1:JG))
p(:,iter) = p(:,iter)/(norm2(Mp(:,iter)))
Mp(:,iter) = Mp(:,iter)/(norm2(Mp(:,iter)))
!Update solution and residual
!call gemv(Mp, r, y [,alpha][,beta] [,trans])
IF(iter.EQ.1) THEN
SUM1 = dot_product(Mp(:,iter),r)
ELSE
DO I=1,MASK
SUM1 = SUM (Mp(:,MASK) * r)
ENDDO
ENDIF
IF ((iter.EQ.1) THEN)
alpha2 = SUM1/SUM(Mp(1:iter,:) * Mp(:,1:iter))
x = x + alpha2 * p(:,iter)
r = r - alpha2 * Mp(:,iter) ! where is the *(Cnew - C)?
!Check convergence
normr=norm2(r)
norm
!fprintf('norm(r) = !g iter = !gNewLine',normr,iter+1);
if (normr < tol) then
!fprintf('GCR Converged, iter = !gNewLine',iter+1);
converged=1
exit
end if
END DO !for loop
if (normr > tol) then
write(*,*) 'GCR SOLUTION DID NOT CONVERGE!'
converged=0
end if
RETURN
END subroutine gcr_Mfree
Вот основные проблемы преобразования формы MATLAB в FORTRAN:
Я не знаю, как найти норму vecto rs и матрицы для условий сходимости? В стандарте 2008 года есть внутреннее заявление (nomr2).
Я не знаю, как правильно обновить код для расчета p, mp, alpha2?
альфа2 уравнение:
alpha2 = r'*Mp(:,iter))/(Mp(iter,:) * Mp(:,iter))
alpha2 definition in MATLAB when iter = 2
alpha2 = (Nx1)'*(Nx2)/(Nx2)' * (Nx2)
(1x2) / (2x2)
(1x2)
Я нашел gemv функции, и это звено, которое описывает вектор для матричных операций массива.
Пожалуйста, более конкретно ... Вектор представляет собой 1D массив, понятие строк и столбцов не применяется здесь! Для 2D-массива вы можете иметь массивы 1xN и Nx1, которые будут соответствовать вектору строки и столбца соответственно. Но тогда вы снова можете использовать 'transpose()'. 'forall' - это просто конструктор/оператор для назначения значений массиву. –
@AlexanderVogt RHS моего уравнения - это скаляр, который является массивом нулевого ранга. LHS состоит из сложного векторного умножения, которое должно давать скаляр в конце. В MATLAB просто я могу использовать ** ** или ** транспонирование **. Но в Фортране это чертовски сложно ... – Vahid
Разве это не только '' dot_product() ') (https://gcc.gnu.org/onlinedocs/gfortran/DOT_005fPRODUCT.html)? –