title: How to properly convert a MATLAB code to Fortran 95 or above?
body: <p>I have written a Matlab code for CGR method which works fine. Here it is: </p>
  
  <p><strong>MATLAB CODE</strong></p>
  
  <pre><code>% 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 &lt; tol)
  %fprintf('GCR Converged, iter = %g\n',iter+1);
  converged=1;
  break;
  end
  
  end %for loop
  
  if(normr &gt; tol)
  fprintf('SOLUTION DID NOT CONVERGE!!');
  converged=0;
  end
  </code></pre>
  
  <p>And here is my last version of Fortran code :</p>
  
  <pre><code>    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  !&lt;
          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 &lt; tol) then
                  !fprintf('GCR Converged, iter = !gNewLine',iter+1);
                  converged=1
                  exit
              end if
  
  
              END DO !for loop
  
              if (normr &gt; tol) then
              write(*,*) 'GCR SOLUTION DID NOT CONVERGE!'
              converged=0
              end if
  
              RETURN
              END subroutine gcr_Mfree
  </code></pre>
  
  <p>Here are the main problems of converting form MATLAB to FORTRAN :</p>
  
  <ol>
  <li><p>I don't know how should I find the norm of vectors and matrices for converging conditions? There is an intrinsic statement (nomr2) at 2008 standard. 
  <a href="https://gcc.gnu.org/onlinedocs/gfortran/NORM2.html" rel="nofollow">norm2 statement in 2008 standard</a></p></li>
  <li><p>I don't know how to appropriately update code to calculate p, mp, alpha2 ?</p></li>
  </ol>
  
  <p>alpha2 equation :</p>
  
  <pre><code> alpha2 = r'*Mp(:,iter)) / (Mp(iter,:) * Mp(:,iter))
  
   alpha2 definition in MATLAB when iter = 2
    alpha2 = (Nx1)'*(Nx2) / (Nx2)' * (Nx2)
                 (1x2)    /     (2x2)
                        (1x2)
  </code></pre>
  
  <p>I found <strong>gemv</strong> function and this link which describes the vector to matrix array operations. </p>
  
  <p><a href="https://software.intel.com/en-us/node/520750" rel="nofollow">gemv function</a> </p>