Program b : gdadcn, gdadgd, gddvgd, gdmpcn, gdmpgd, gdsqrt, sqadsq,

      program progb
      implicit none

      real rms
      external rms

      integer i,j
      real matax(3,4),matay(3,4),matbx(3,4),matby(3,4)
      real matc(3,4),resultx(3,4),resulty(3,4),length(3,4)
      real result2x(3,4),result2y(3,4)
# initialize mata, matb and matc
      do i=1,3
         do j=1,4
            matax(i,j)=float(i-3.11*j)
            matay(i,j)=float(i+2.1*j)
            matbx(i,j)=float(3.45*i*j)
            matby(i,j)=float(i)/float(j)
            matc(i,j)=float(-i+2.13*j)
         enddo
      enddo

      print *,'The calculation that will be done is:'
      print *,'Result=(5*A+matrix of (1.0,3.0)+B) *** C '
      print *,'Result=Result/Length(Result)*10'
      print *,'A,B and Result are vector matrices'
      print *,'C and Length are scalar matrices'
      print *,'*** means element-wise multiplication by a scalar'
      print *,''
      print *,'matrix A equals:'
      do j=1,4
         print *,'(',matax(1,j),',',matay(1,j),' )  (',
     *    matax(2,j),',',matay(2,j),' )  (',
     *    matax(3,j),',',matay(3,j),' )' 
      enddo
      print *,''
      print *,'matrix B equals:'
      do j=1,4
         print *,'(',matbx(1,j),',',matby(1,j),' )  (',
     *    matbx(2,j),',',matby(2,j),' )  (',
     *    matbx(3,j),',',matby(3,j),' )' 
      enddo   
      print *,''
 
      print *,'matrix C equals:'
      do j=1,4
         print *,matc(1,j),' ',matc(2,j),' ',matc(3,j)
      enddo
      print *,''

# make the calculation with two method

# compute 5*A+matrix of (1.0,3.0)
      call gdmpcn(resultx,matax,5.0,3,4,0)
      call gdmpcn(resulty,matay,5.0,3,4,0)
      call gdadcn(resultx,resultx,1.0,3,4,0)
      call gdadcn(resulty,resulty,3.0,3,4,0)

# compute (Result+B) *** C
      call gdadgd(resultx,resultx,matbx,1.0,1.0,3,4,0)
      call gdadgd(resulty,resulty,matby,1.0,1.0,3,4,0)
      call gdmpgd(resultx,resultx,matc,1.0,3,4,0)
      call gdmpgd(resulty,resulty,matc,1.0,3,4,0)
      
# compute Length(Result), a 3x4 matrix     
      call sqadsq(length,resultx,resulty,1.0,1.0,3,4,0)
      call gdsqrt(length,length,1.0,3,4,0)

# normalize vectors of Result and multiply them by 10
      call gddvgd(resultx,resultx,length,10.0,3,4,0)
      call gddvgd(resulty,resulty,length,10.0,3,4,0)

# print Result
      print *,'matrix Result equals:'
      do j=1,4
         print *,'(',resultx(1,j),',',resulty(1,j),' )  (',
     *    resultx(2,j),',',resulty(2,j),' )  (',
     *    resultx(3,j),',',resulty(3,j),' )' 
      enddo   
      print *,''


# second method

      do i=1,3
         do j=1,4
            result2x(i,j)=(5*matax(i,j)+1.0+matbx(i,j))*matc(i,j)
            result2y(i,j)=(5*matay(i,j)+3.0+matby(i,j))*matc(i,j)
            length(i,j)=sqrt(result2x(i,j)**2+result2y(i,j)**2)
            if (length(i,j).ne.0) then
               result2x(i,j)=10*result2x(i,j)/length(i,j)
               result2y(i,j)=10*result2y(i,j)/length(i,j)
            endif   
         enddo
      enddo   
      print *,'matrix Result2 equals:'
      do j=1,4
         print *,'(',result2x(1,j),',',result2y(1,j),' )  (',
     *    result2x(2,j),',',result2y(2,j),' )  (',
     *    result2x(3,j),',',result2y(3,j),' )' 
      enddo  

      print *,'The rms difference between the two methods is:',
     * rms(resultx,result2x,3,4)+rms(resulty,result2y,3,4) 

      stop
      end

# computes the rms of the difference of arrays a and b
      real function rms(a,b,n,m)
      real a(n,m),b(n,m),sum
      integer i,j,n,m
      sum=0
      do i=1,n
         do j=1,m
            sum=sum+(a(i,j)-b(i,j))**2
         enddo
      enddo
      rms=sqrt(sum/n/m)
      return
      end