!
!                                **********************************************
!                                *             MODULE peuc                    *
!                                *  R. J. Purser, NOAA/NCEP/EMC      Oct 2005 * 
!                                *                              18th May 2012 *
!                                *  jim.purser@noaa.gov                       *
!                                *                                            *
!                                **********************************************
!
!  Euclidean geometry, geometric (stereographic) projections,
!              related transformations (Mobius).
! Package for handy vector and matrix operations in Euclidean geometry.
! This package is primarily intended for 3D operations and three of the
! functions (Cross_product, Triple_product and Axial) do not possess simple
! generalizations to a generic number N of dimensions. The others, while
! admitting such N-dimensional generalizations, have not all been provided
! with such generic forms here at the time of writing, though some of these 
! may be added at a future date.
!
!   FUNCTION:
! Normalized:     Normalized version of given real vector
! Cross_product:  Vector cross-product of the given 2 vectors
! Outer_product:  outer-product matrix of the given 2 vectors
! Triple_product: Scalar triple product of given 3 vectors
! Det:            Determinant of given matrix
! Axial:          Convert axial-vector <--> 2-form (antisymmetric matrix)
! Diag:           Diagnl of given matrix, or diagonal matrix of given elements
! Trace:          Trace of given matrix
! Identity:       Identity 3*3 matrix, or identity n*n matrix for a given n
! Sarea:          Spherical area subtended by three vectors
! Huarea:         Spherical area subtended by right-angled spherical triangle
!   SUBROUTINE:
! Gram:           Right-handed orthogonal basis and rank, nrank. The first 
!                 nrank basis vectors span the column range of matrix given,
!                 OR  ("plain" version) simple unpivoted Gram-Schmidt of a
!                 square matrix.
!
! In addition, we include routines that relate to stereographic projections
! and some associated mobius transformation utilities, since these complex
! operations have a strong geomoetrical flavor.
!
! DIRECT DEPENDENCIES
! Libraries[their Modules]: pmat[pmat]
! Additional Modules      : pkind
!
!============================================================================
module peuc
!============================================================================
use pkind, only: sp,dp,dpc
implicit none
private
public:: normalized,cross_product,outer_product,triple_product,det,axial, &
         diag,trace,identity,sarea,huarea,gram,rowops,corral, &
         ctoz,ztoc,setmobius, &
         mobius,mobiusi

interface normalized;module procedure normalized_s,normalized_d;end interface
interface cross_product
   module procedure cross_product_s,cross_product_d;            end interface
interface outer_product
   module procedure outer_product_s,outer_product_d,outer_product_i
                                                                end interface
interface triple_product
   module procedure triple_product_s,triple_product_d;          end interface
interface det; module procedure det_s,det_d,det_i;              end interface
interface axial
   module procedure axial3_s,axial3_d,axial33_s,axial33_d;      end interface
interface diag
   module procedure diagn_s,diagn_d,diagn_i,diagnn_s,diagnn_d,diagnn_i
                                                                end interface
interface trace;    module procedure trace_s,trace_d,trace_i;   end interface
interface identity; module procedure identity_i,identity3_i;    end interface
interface huarea;   module procedure huarea_s,huarea_d;         end interface
interface sarea;    module procedure sarea_s,sarea_d;           end interface
interface gram
   module procedure gram_s,gram_d,plaingram_s,plaingram_d,rowgram
                                                                end interface
interface rowops;   module procedure rowops;                    end interface
interface corral;   module procedure corral;                    end interface
interface ctoz;     module procedure ctoz;                      end interface
interface ztoc;     module procedure ztoc,ztocd;                end interface
interface setmobius;module procedure setmobius,zsetmobius;      end interface
interface mobius;   module procedure zmobius,cmobius;           end interface
interface mobiusi;  module procedure zmobiusi;                  end interface

contains

!=============================================================================
function normalized_s(a)result(b)!                                [normalized]
!=============================================================================
real(sp),dimension(:),intent(IN):: a
real(sp),dimension(size(a))     :: b
real(sp)                        :: s
s=sqrt(dot_product(a,a)); if(s==0)then; b=0;else;b=a/s;endif
end function normalized_s
!=============================================================================
function normalized_d(a)result(b)!                                [normalized]

!=============================================================================
real(dp),dimension(:),intent(IN):: a
real(dp),dimension(size(a))     :: b
real(dp)                        :: s
s=sqrt(dot_product(a,a)); if(s==0)then; b=0;else;b=a/s;endif
end function normalized_d

!=============================================================================
function cross_product_s(a,b)!                                 [cross_product]
!=============================================================================
real(sp),dimension(3)           :: cross_product_s
real(sp),dimension(3),intent(in):: a,b
cross_product_s(1)=a(2)*b(3)-a(3)*b(2)
cross_product_s(2)=a(3)*b(1)-a(1)*b(3)
cross_product_s(3)=a(1)*b(2)-a(2)*b(1)
end function cross_product_s
!=============================================================================
function cross_product_d(a,b)result(cross_product)!            [cross_product]
!=============================================================================
real(dp),dimension(3)           :: cross_product
real(dp),dimension(3),intent(in):: a,b
cross_product(1)=a(2)*b(3)-a(3)*b(2)
cross_product(2)=a(3)*b(1)-a(1)*b(3)
cross_product(3)=a(1)*b(2)-a(2)*b(1)
end function cross_product_d

!=============================================================================
function outer_product_s(a,b)result(c)!                        [outer_product]
!=============================================================================
real(sp),dimension(:),  intent(in ):: a
real(sp),dimension(:),  intent(in ):: b
real(sp),DIMENSION(size(a),size(b)):: c
integer                            :: nb,i
nb=size(b)
do i=1,nb; c(:,i)=a*b(i); enddo
end function outer_product_s
!=============================================================================
function outer_product_d(a,b)result(c)!                        [outer_product]
!=============================================================================
real(dp),dimension(:),  intent(in ):: a
real(dp),dimension(:),  intent(in ):: b
real(dp),dimension(size(a),size(b)):: c
integer                            :: nb,i
nb=size(b)
do i=1,nb; c(:,i)=a*b(i); enddo
end function outer_product_d
!=============================================================================
function outer_product_i(a,b)result(c)!                        [outer_product]
!=============================================================================
integer,dimension(:),  intent(in ):: a
integer,dimension(:),  intent(in ):: b
integer,dimension(size(a),size(b)):: c
integer                            :: nb,i
nb=size(b)
do i=1,nb; c(:,i)=a*b(i); enddo
end function outer_product_i

!=============================================================================
function triple_product_s(a,b,c)result(tripleproduct)!        [triple_product]
!=============================================================================
real(sp),dimension(3),intent(IN ):: a,b,c
real(sp)                         :: tripleproduct
tripleproduct=dot_product( cross_product(a,b),c )
end function triple_product_s
!=============================================================================
function triple_product_d(a,b,c)result(tripleproduct)!        [triple_product]
!=============================================================================
real(dp),dimension(3),intent(IN ):: a,b,c
real(dp)                         :: tripleproduct
tripleproduct=dot_product( cross_product(a,b),c )
end function triple_product_d

!=============================================================================
function det_s(a)result(det)!                                            [det]
!=============================================================================
real(sp),dimension(:,:),intent(IN )    ::a
real(sp)                               :: det
real(sp),dimension(size(a,1),size(a,1)):: b
integer                                :: n,nrank
n=size(a,1)
if(n==3)then
   det=triple_product(a(:,1),a(:,2),a(:,3))
else
   call gram(a,b,nrank,det)
   if(nrank<n)det=0
endif
end function det_s
!=============================================================================
function det_d(a)result(det)!                                            [det]
!=============================================================================
real(dp),dimension(:,:),intent(IN )    ::a
real(dp)                               :: det
real(dp),dimension(size(a,1),size(a,1)):: b
integer                                :: n,nrank
n=size(a,1)
if(n==3)then
   det=triple_product(a(:,1),a(:,2),a(:,3))
else
   call gram(a,b,nrank,det)
   if(nrank<n)det=0
endif
end function det_d
!=============================================================================
function det_i(a)result(idet)!                                           [det]
!=============================================================================
integer, dimension(:,:),intent(IN )    :: a
integer                                :: idet
real(dp),dimension(size(a,1),size(a,2)):: b
real(dp)                               :: bdet
b=a; bdet=det(b); idet=nint(bdet)
end function det_i

!=============================================================================
function axial3_s(a)result(b)!                                         [axial]
!=============================================================================
real(sp),dimension(3),intent(IN ):: a
real(sp),dimension(3,3)          :: b
b=0;b(3,2)=a(1);b(1,3)=a(2);b(2,1)=a(3);b(2,3)=-a(1);b(3,1)=-a(2);b(1,2)=-a(3)
end function axial3_s
!=============================================================================
function axial3_d(a)result(b)!                                         [axial]
!=============================================================================
real(dp),dimension(3),intent(IN ):: a
real(dp),dimension(3,3)          :: b
b=0;b(3,2)=a(1);b(1,3)=a(2);b(2,1)=a(3);b(2,3)=-a(1);b(3,1)=-a(2);b(1,2)=-a(3)
end function axial3_d
!=============================================================================
function axial33_s(b)result(a)!                                        [axial]
!=============================================================================
real(sp),dimension(3,3),intent(IN ):: b
real(sp),dimension(3)              :: a
a(1)=(b(3,2)-b(2,3))/2; a(2)=(b(1,3)-b(3,1))/2; a(3)=(b(2,1)-b(1,2))/2
end function axial33_s
!=============================================================================
function axial33_d(b)result(a)!                                        [axial]
!=============================================================================
real(dp),dimension(3,3),intent(IN ):: b
real(dp),dimension(3)              :: a
a(1)=(b(3,2)-b(2,3))/2; a(2)=(b(1,3)-b(3,1))/2; a(3)=(b(2,1)-b(1,2))/2
end function axial33_d

!=============================================================================
function diagn_s(a)result(b)!                                           [diag]
!=============================================================================
real(sp),dimension(:),intent(IN )  :: a
real(sp),dimension(size(a),size(a)):: b
integer                            :: n,i
n=size(a)
b=0; do i=1,n; b(i,i)=a(i); enddo
end function diagn_s
!=============================================================================
function diagn_d(a)result(b)!                                           [diag]
!=============================================================================
real(dp),dimension(:),intent(IN )  :: a
real(dp),dimension(size(a),size(a)):: b
integer                            :: n,i
n=size(a)
b=0; do i=1,n; b(i,i)=a(i); enddo
end function diagn_d
!=============================================================================
function diagn_i(a)result(b)!                                           [diag]
!=============================================================================
integer,dimension(:),intent(IN )  :: a
integer,dimension(size(a),size(a)):: b
integer                            :: n,i
n=size(a)
b=0; do i=1,n; b(i,i)=a(i); enddo
end function diagn_i
!=============================================================================
function diagnn_s(b)result(a)!                                          [diag]
!=============================================================================
real(sp),dimension(:,:),intent(IN ):: b
real(sp),dimension(size(b,1))      :: a
integer                            :: n,i
n=size(b,1)
do i=1,n; a(i)=b(i,i); enddo
end function diagnn_s
!=============================================================================
function diagnn_d(b)result(a)!                                          [diag]
!=============================================================================
real(dp),dimension(:,:),intent(IN ):: b
real(dp),dimension(size(b,1))      :: a
integer                            :: n,i
n=size(b,1)
do i=1,n; a(i)=b(i,i); enddo
end function diagnn_d
!=============================================================================
function diagnn_i(b)result(a)!                                          [diag]
!=============================================================================
integer,dimension(:,:),intent(IN ):: b
integer,dimension(size(b,1))      :: a
integer                            :: n,i
n=size(b,1)
do i=1,n; a(i)=b(i,i); enddo
end function diagnn_i

!=============================================================================
function trace_s(b)result(s)!                                          [trace]
!=============================================================================
real(sp),dimension(:,:),intent(IN ):: b
real(sp)                           :: s
s=sum(diag(b))
end function trace_s
!=============================================================================
function trace_d(b)result(s)!                                          [trace]
!=============================================================================
real(dp),dimension(:,:),intent(IN ):: b
real(dp)                           :: s
s=sum(diag(b))
end function trace_d
!=============================================================================
function trace_i(b)result(s)!                                         [trace]
!=============================================================================
integer,dimension(:,:),intent(IN ):: b
integer                           :: s
s=sum(diag(b))
end function trace_i

!=============================================================================
function identity_i(n)result(a)!                                    [identity]
!=============================================================================
integer,intent(IN )   :: n
integer,dimension(n,n):: a
integer               :: i
a=0; do i=1,n; a(i,i)=1; enddo
end function identity_i
!=============================================================================
function identity3_i()result(a)!                                    [identity]
!=============================================================================
integer,dimension(3,3):: a
integer               :: i
a=0; do i=1,3; a(i,i)=1; enddo
end function identity3_i

!=============================================================================
function huarea_s(sa,sb)result(area)!                                 [huarea]
!=============================================================================
implicit none
real(sp),intent(IN ):: sa,sb
real(sp)            :: area
real(sp)            :: ca,cb
!-----------------------------------------------------------------------------
ca=sqrt(1-sa**2)
cb=sqrt(1-sb**2)
area=asin(sa*sb/(1+ca*cb))
end function huarea_s

!=============================================================================
function huarea_d(sa,sb)result(area)!                                 [huarea]
!=============================================================================
implicit none
real(dp),intent(IN ):: sa,sb
real(dp)            :: area
real(dp)            :: ca,cb
!-----------------------------------------------------------------------------
ca=sqrt(1-sa**2)
cb=sqrt(1-sb**2)
area=asin(sa*sb/(1+ca*cb))
end function huarea_d


!=============================================================================
function sarea_s(v1,v2,v3)result(area)!                                [sarea]
!=============================================================================
! Compute the area of the spherical triangle, {v1,v2,v3}, measured in the
! right-handed sense, by dropping a perpendicular to u0 on the longest side so
! that the area becomes the sum of areas of the two simpler right-angled
! triangles.
!=============================================================================
real(sp),dimension(3),intent(IN ):: v1,v2,v3
real(sp)                         :: area
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
real(sp)                         :: s123,a1,a2,b,d1,d2,d3
real(sp),dimension(3)            :: u0,u1,u2,u3,x,y
!=============================================================================
area=0
u1=normalized(v1); u2=normalized(v2); u3=normalized(v3)
s123=triple_product(u1,u2,u3)
if(s123==0)return

d1=dot_product(u3-u2,u3-u2)
d2=dot_product(u1-u3,u1-u3)
d3=dot_product(u2-u1,u2-u1)

! Triangle that is not degenerate. Cyclically permute, so side 3 is longest:
if(d3<d1 .or. d3<d2)call cyclic(u1,u2,u3,d1,d2,d3)
if(d3<d1 .or. d3<d2)call cyclic(u1,u2,u3,d1,d2,d3)
y=normalized( cross_product(u1,u2) )
b=dot_product(y,u3)
u0=normalized( u3-y*b )
x=cross_product(y,u0) 
a1=-dot_product(x,u1-u0); a2= dot_product(x,u2-u0)
area=huarea(a1,b)+huarea(a2,b)

contains
!-----------------------------------------------------------------------------
   subroutine cyclic(u1,u2,u3,d1,d2,d3)
!-----------------------------------------------------------------------------
   real(sp),dimension(3),intent(INOUT):: u1,u2,u3
   real(sp),             intent(INOUT):: d1,d2,d3
   real(sp),dimension(3)              :: ut
   real(sp)                           :: dt
   dt=d1; d1=d2; d2=d3; d3=dt
   ut=u1; u1=u2; u2=u3; u3=ut
   end subroutine cyclic
end function sarea_s

!=============================================================================
function sarea_d(v1,v2,v3)result(area)!                                [sarea]
!=============================================================================
real(dp),dimension(3),intent(IN ):: v1,v2,v3
real(dp)                         :: area
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
real(dp)                         :: s123,a1,a2,b,d1,d2,d3
real(dp),dimension(3)            :: u0,u1,u2,u3,x,y
!=============================================================================
area=0
u1=normalized(v1); u2=normalized(v2); u3=normalized(v3)
s123=triple_product(u1,u2,u3)
if(s123==0)return

d1=dot_product(u3-u2,u3-u2)
d2=dot_product(u1-u3,u1-u3)
d3=dot_product(u2-u1,u2-u1)

! Triangle that is not degenerate. Cyclically permute, so side 3 is longest:
if(d3<d1 .or. d3<d2)call cyclic(u1,u2,u3,d1,d2,d3)
if(d3<d1 .or. d3<d2)call cyclic(u1,u2,u3,d1,d2,d3)
y=normalized( cross_product(u1,u2) )
b=dot_product(y,u3)
u0=normalized( u3-y*b )
x=cross_product(y,u0) 
a1=-dot_product(x,u1-u0); a2= dot_product(x,u2-u0)
area=huarea(a1,b)+huarea(a2,b)

contains
!-----------------------------------------------------------------------------
   subroutine cyclic(u1,u2,u3,d1,d2,d3)
!-----------------------------------------------------------------------------
   real(dp),dimension(3),intent(INOUT):: u1,u2,u3
   real(dp),             intent(INOUT):: d1,d2,d3
   real(dp),dimension(3)              :: ut
   real(dp)                           :: dt
   dt=d1; d1=d2; d2=d3; d3=dt
   ut=u1; u1=u2; u2=u3; u3=ut
   end subroutine cyclic
end function sarea_d

!=============================================================================
subroutine gram_s(as,b,nrank,det)!                                      [gram]
!=============================================================================
real(sp),dimension(:,:),intent(IN )      :: as
real(sp),dimension(:,:),intent(OUT)      :: b
integer,                intent(OUT)      :: nrank
real(sp),               intent(OUT)      :: det
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
real(sp),parameter                       :: crit=1.e-5_sp
real(sp),dimension(size(as,1),size(as,2)):: a
real(sp),dimension(size(as,2),size(as,1)):: ab
real(sp),dimension(size(as,1))           :: tv,w
real(sp)                                 :: val,s,vcrit
integer                                  :: i,j,k,l,m,n
integer,dimension(2)                     :: ii
!=============================================================================
n=size(as,1)
m=size(as,2)
if(n/=size(b,1) .or. n/=size(b,2))stop 'In gram; incompatible dimensions'
a=as
b=identity(n)
det=1
val=maxval(abs(a))
if(val==0)then
   nrank=0
   return
endif
vcrit=val*crit
nrank=min(n,m)
do k=1,n
   if(k>nrank)exit
   ab(k:m,k:n)=matmul( transpose(a(:,k:m)),b(:,k:n) )
   ii =maxloc( abs( ab(k:m,k:n)) )+k-1
   val=maxval( abs( ab(k:m,k:n)) )
   if(val<=vcrit)then
      nrank=k-1
      exit
   endif
   i=ii(1)
   j=ii(2)
   tv=b(:,j)
   b(:,j)=-b(:,k)
   b(:,k)=tv
   tv=a(:,i)
   a(:,i)=-a(:,k)
   a(:,k)=tv
   w(k:n)=matmul( transpose(b(:,k:n)),tv )
   b(:,k)=matmul(b(:,k:n),w(k:n) )
   s=dot_product(b(:,k),b(:,k))
   s=sqrt(s)
   if(w(k)<0)s=-s
   det=det*s
   b(:,k)=b(:,k)/s
   do l=k,n
      do j=l+1,n
         s=dot_product(b(:,l),b(:,j))
         b(:,j)=normalized( b(:,j)-b(:,l)*s )
      enddo
   enddo
enddo
end subroutine gram_s
   
!=============================================================================
subroutine gram_d(as,b,nrank,det)!                                      [gram]
!=============================================================================
real(dp),dimension(:,:),intent(IN )      :: as
real(dp),dimension(:,:),intent(OUT)      :: b
integer,                intent(OUT)      :: nrank
real(dp),               intent(OUT)      :: det
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
real(dp),parameter                       :: crit=1.e-9_dp
real(dp),dimension(size(as,1),size(as,2)):: a
real(dp),dimension(size(as,2),size(as,1)):: ab
real(dp),dimension(size(as,1))           :: tv,w
real(dp)                                 :: val,s,vcrit
integer                                  :: i,j,k,l,m,n
integer,dimension(2)                     :: ii
!=============================================================================
n=size(as,1)
m=size(as,2)
if(n/=size(b,1) .or. n/=size(b,2))stop 'In gram; incompatible dimensions'
a=as
b=identity(n)
det=1
val=maxval(abs(a))
if(val==0)then
   nrank=0
   return
endif
vcrit=val*crit
nrank=min(n,m)
do k=1,n
   if(k>nrank)exit
   ab(k:m,k:n)=matmul( transpose(a(:,k:m)),b(:,k:n) )
   ii =maxloc( abs( ab(k:m,k:n)) )+k-1
   val=maxval( abs( ab(k:m,k:n)) )
   if(val<=vcrit)then
      nrank=k-1
      exit
   endif
   i=ii(1)
   j=ii(2)
   tv=b(:,j)
   b(:,j)=-b(:,k)
   b(:,k)=tv
   tv=a(:,i)
   a(:,i)=-a(:,k)
   a(:,k)=tv
   w(k:n)=matmul( transpose(b(:,k:n)),tv )
   b(:,k)=matmul(b(:,k:n),w(k:n) )
   s=dot_product(b(:,k),b(:,k))
   s=sqrt(s)
   if(w(k)<0)s=-s
   det=det*s
   b(:,k)=b(:,k)/s
   do l=k,n
      do j=l+1,n
         s=dot_product(b(:,l),b(:,j))
         b(:,j)=normalized( b(:,j)-b(:,l)*s )
      enddo
   enddo
enddo
end subroutine gram_d
   
!=============================================================================
subroutine plaingram_s(b,nrank)!                                        [gram]
!=============================================================================
! A "plain" (unpivoted) version of Gram-Schmidt, for square matrices only.
real(sp),dimension(:,:),intent(INOUT)    :: b
integer,                intent(  OUT)    :: nrank
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
real(sp),parameter                       :: crit=1.e-5_sp
real(sp)                                 :: val,vcrit
integer                                  :: j,k,n
!=============================================================================
n=size(b,1); if(n/=size(b,2))stop 'In gram; matrix needs to be square'
val=maxval(abs(b))
nrank=0
if(val==0)then
   b=0
   return
endif
vcrit=val*crit
do k=1,n
   val=sqrt(dot_product(b(:,k),b(:,k)))
   if(val<=vcrit)then
      b(:,k:n)=0
      return
   endif
   b(:,k)=b(:,k)/val
   nrank=k
   do j=k+1,n
      b(:,j)=b(:,j)-b(:,k)*dot_product(b(:,k),b(:,j))
   enddo
enddo
end subroutine plaingram_s

!=============================================================================
subroutine plaingram_d(b,nrank)!                                        [gram]
!=============================================================================
! A "plain" (unpivoted) version of Gram-Schmidt, for square matrices only.
real(dp),dimension(:,:),intent(INOUT)    :: b
integer,                intent(  OUT)    :: nrank
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
real(dp),parameter                       :: crit=1.e-9_dp
real(dp)                                 :: val,vcrit
integer                                  :: j,k,n
!=============================================================================
n=size(b,1); if(n/=size(b,2))stop 'In gram; matrix needs to be square'
val=maxval(abs(b))
nrank=0
if(val==0)then
   b=0
   return
endif
vcrit=val*crit
do k=1,n
   val=sqrt(dot_product(b(:,k),b(:,k)))
   if(val<=vcrit)then
      b(:,k:n)=0
      return
   endif
   b(:,k)=b(:,k)/val
   nrank=k
   do j=k+1,n
      b(:,j)=b(:,j)-b(:,k)*dot_product(b(:,k),b(:,j))
   enddo
enddo
end subroutine plaingram_d

!=============================================================================
subroutine rowgram(m,n,a,ipiv,tt,b,rank)!                               [gram]
!=============================================================================
! Without changing (tall) rectangular input matrix a, perform pivoted gram-
! schmidt operations to orthogonalize the rows, until rows that remain become
! negligible. Record the pivoting sequence in ipiv, and the row-normalization
! in tt(j,j) and the row-orthogonalization in tt(i,j), for i>j. Note that
! tt(i,j)=0 for i<j (tt is truncated lower triangular). The orthonormalized
! rows are returned in square array b, which is complete eben when the
! effective rank < n. 
! The recorded row operations can be repeated on independent column vectors
! through the use of subroutine ROWOPS (in this module).
! It is recommended to rescale the original matrix A via a call to CORRAL
! (in this module) because the negligibility criterion depends upon an
! "epsilon" value that is fixed (10**(-13)) and assumes elements of a are
! never too different in magnitude from unity, unless they are actually zero.
!=============================================================================
use pkind
implicit none
integer,                intent(IN ):: m,n
real(dp),dimension(m,n),intent(in ):: a
integer, dimension(n),  intent(out):: ipiv
real(dp),dimension(m,n),intent(out):: tt
real(dp),dimension(n,n),intent(out):: b
integer,                intent(out):: rank
!-----------------------------------------------------------------------------
real(dp),parameter      :: eps=1.e-13_dp,epss=eps**2
real(dp),dimension(m,n) :: aa
real(dp),dimension(n)   :: rowv
real(dp),dimension(m)   :: p
real(dp)                :: maxp,nepss
integer,dimension(1)    :: jloc
integer                 :: i,ii,iii,j,maxi
!=============================================================================
if(m<n)stop 'In rowgram; this routines needs m>=n please'
nepss=n*epss
rank=n
aa=a
tt=0
do ii=1,n

! At this stage, all rows less than ii are already orthonormalized and are
! orthogonal to all rows at and beyond ii. Find the norms of these lower
! rows and pivot the largest of them into position ii:
   maxp=0
   maxi=ii
   do i=ii,m
      p(i)=dot_product(aa(i,:),aa(i,:))
      if(p(i)>maxp)then
         maxp=p(i)
         maxi=i
      endif
   enddo
   if(maxp<nepss)then !<- End of gram process; clean up and return
      b=0
      b(1:ii-1,:)=aa(1:ii-1,:)
! fill the remaining rows, ii:n, of b with remaining orthonormal rows at random
      do iii=ii,n
! find the column of b for which the maximum element is the smallest:
         do j=1,n
            rowv(j)=maxval(abs(b(1:iii-1,j)))
         enddo
         jloc=minloc(rowv)
         j=jloc(1)
         b(iii,j)=1
         do i=1,iii-1
            maxp=dot_product(b(i,:),b(iii,:))
            b(iii,:)=b(iii,:)-b(i,:)*maxp
         enddo
         maxp=sqrt(dot_product(b(iii,:),b(iii,:)))
         b(iii,:)=b(iii,:)/maxp
      enddo
      rank=ii-1
      return
   endif
   
   ipiv(ii)=maxi
   if(maxi/=ii)then
      rowv      =aa(ii,  :)
      aa(ii,  :)=aa(maxi,:)
      aa(maxi,:)=rowv
   endif
   maxp=sqrt(maxp)
   tt(ii,ii)=maxp
   aa(ii,:)=aa(ii,:)/maxp
! Adjust all rows below to make them orthogonal to new row ii  
   do i=ii+1,m
      maxp=dot_product(aa(ii,:),aa(i,:))
      tt(i,ii)=maxp
      aa(i,:)=aa(i,:)-aa(ii,:)*maxp
   enddo
enddo
b=aa(1:n,:)
end subroutine rowgram

!=============================================================================
subroutine rowops(m,n,ipiv,tt,v,vv)!                                  [rowops]
!=============================================================================
! Apply the row-operations, implied by ipiv and tt returned by rowgram, to
! the single column vector, v, to produce the transformed vector vv.
!=============================================================================
use pkind
implicit none
integer,                intent(in ):: m,n
integer, dimension(n),  intent(in ):: ipiv
real(dp),dimension(m,n),intent(in ):: tt
real(dp),dimension(m),  intent(in ):: v
real(dp),dimension(m),  intent(out):: vv
!-----------------------------------------------------------------------------
integer  :: i,j,k
real(dp) :: p
!=============================================================================
vv=v
do j=1,n
   k=ipiv(j)
   if(k/=j)then
      p=vv(j)
      vv(j)=vv(k)
      vv(k)=p
   endif
   vv(j)=vv(j)/tt(j,j)
   do i=j+1,m
      vv(i)=vv(i)-vv(j)*tt(i,j)
   enddo
enddo
end subroutine rowops
   
!=============================================================================
subroutine corral(m,n,mask,a,d,aa,e)!                                 [corral]
!=============================================================================
! Find positive diagonals D and E and a Lagrange multiplier F that minimize
! the row-sum +column-sum of masked terms, 
! (D_i +log(|A_ij|) +E_j)^2
! subject to the single constraint, sum_j E_j =0, where the mask permits
! only nonnegligible A_ij to participate in the quadratic quantities. 
! Once a solution for D and E is found, return their exponentials, d and e,
! together with the rescaled matrix aa such that a = d.aa.e when d and e are
! interpretted as diagonal matrices.
!=============================================================================
use pkind
use pmat
implicit none
integer,                intent(in ):: m,n
logical, dimension(m,n),intent(in ):: mask
real(dp),dimension(m,n),intent(in ):: a
real(dp),dimension(m  ),intent(out):: d
real(dp),dimension(m,n),intent(out):: aa
real(dp),dimension(  n),intent(out):: e
!-----------------------------------------------------------------------------
real(dp),dimension(0:m+n,0:m+n):: g
real(dp),dimension(0:m+n)      :: h
integer                        :: i,j,k,nh
!=============================================================================
nh=1+m+n
aa=0
do j=1,n
do i=1,m
   if(mask(i,j))aa(i,j)=log(abs(a(i,j)))
enddo
enddo

h=0
g=0

! Equations on row 0 enforcing Lagrange multiplier F-constraint:
do j=1,n
   k=m+j
   g(0,k)=1
enddo

! Equations on rows 1:m minimizing row sums of quadratic terms:
do i=1,m
   do j=1,n
      k=m+j
      if(mask(i,j))then
         g(i,i)=g(i,i)-1
         g(i,k)=-1
         h(i)=h(i)-aa(i,j)
      endif
   enddo
enddo

! Equations on rows m+1:m+n minimizing col sums subject to constraint
do j=1,n
   k=m+j
   g(k,0)=1
   do i=1,m
      if(mask(i,j))then
         g(k,k)=g(k,k)-1
         g(k,i)=-1
         h(k)=h(k)-aa(i,j)
      endif
   enddo
enddo

! Invert the normal equations:
call inv(g,h)

! Exponentiate the parts that become final scaling diagnonal matrices d and e:
do i=1,m
   d(i)=exp(h(i))
enddo
do j=1,n
   k=m+j
   e(j)=exp(h(k))
enddo

! Compute the rescaled matrix directly:
do j=1,n
do i=1,m
   aa(i,j)=a(i,j)/(d(i)*e(j))
enddo
enddo
end subroutine corral
   
!=============================================================================
! Utility code for various Mobius transformations. If aa1,bb1,cc1,dd1 are
! the coefficients for one transformation, and aa2,bb2,cc2,dd2 are the 
! coefficients for a second one, then the coefficients for the mapping
! of a test point, zz, by aa1 etc to zw, followed by a mapping of zw, by
! aa2 etc to zv, is equivalent to a single mapping zz-->zv by the transformatn
! with coefficients aa3,bb3,cc3,dd3, such that, as 2*2 complex matrices:
!
! [ aa3, bb3 ]   [ aa2, bb2 ]   [ aa1, bb1 ]
! [          ] = [          ] * [          ]
! [ cc3, dd3 ]   [ cc2, dd2 ]   [ cc1, dd1 ] .
!
!  Note that the determinant of these matrices is always +1
!
!=============================================================================
subroutine ctoz(v, z,infz)!                                             [ctoz]
!=============================================================================
real(dp),dimension(3),intent(IN ):: v
complex(dpc),         intent(OUT):: z
logical,              intent(OUT):: infz
!-----------------------------------------------------------------------------
real(dp),parameter:: zero=0,one=1
real(dp)          :: rr,zzpi
!=============================================================================
infz=.false.
z=cmplx(v(1),v(2),dpc)
if(v(3)>0)then
   zzpi=one/(one+v(3))
else
   rr=v(1)**2+v(2)**2
   infz=(rr==zero); if(infz)return ! <- The point is mapped to infinity (90S)
   zzpi=(one-v(3))/rr
endif
z=z*zzpi
end subroutine ctoz
   
!=============================================================================
subroutine ztoc(z,infz, v)!                                             [ztoc]
!=============================================================================
complex(dpc),         intent(IN ):: z
logical,              intent(IN ):: infz
real(dp),dimension(3),intent(OUT):: v
!-----------------------------------------------------------------------------
real(dp),parameter:: zero=0,one=1
real(dp)          :: r,q,rs,rsc,rsbi
!=============================================================================
if(infz)then; v=(/zero,zero,-one/); return; endif
r=real(z); q=aimag(z); rs=r*r+q*q
rsc=one-rs
rsbi=one/(one+rs)
v(1)=2*rsbi*r
v(2)=2*rsbi*q
v(3)=rsc*rsbi
end subroutine ztoc

!=============================================================================
subroutine ztocd(z,infz, v,vd)!                                         [ztoc]
!=============================================================================
! The convention adopted for the complex derivative is that, for a complex
! infinitesimal map displacement, delta_z, the corresponding infinitesimal
! change of cartesian vector position is delta_v given by:
! delta_v = Real(vd*delta_z).
! Thus, by a kind of Cauchy-Riemann relation, Imag(vd)=v CROSS Real(vd).
! THE DERIVATIVE FOR THE IDEAL POINT AT INFINITY HAS NOT BEEN CODED YET!!!
!=============================================================================
complex(dpc),             intent(IN ):: z
logical,                  intent(IN ):: infz
real(dp),dimension(3),    intent(OUT):: v
complex(dpc),dimension(3),intent(OUT):: vd
!-----------------------------------------------------------------------------
real(dp),parameter   :: zero=0,one=1
real(dp)             :: r,q,rs,rsc,rsbi,rsbis
real(dp),dimension(3):: u1,u2
integer              :: i
!=============================================================================
if(infz)then; v=(/zero,zero,-one/); return; endif
r=real(z); q=aimag(z); rs=r*r+q*q
rsc=one-rs
rsbi=one/(one+rs)
rsbis=rsbi**2
v(1)=2*rsbi*r
v(2)=2*rsbi*q
v(3)=rsc*rsbi
u1(1)=2*(one+q*q-r*r)*rsbis
u1(2)=-4*r*q*rsbis
u1(3)=-4*r*rsbis
u2=cross_product(v,u1)
do i=1,3
   vd(i)=cmplx(u1(i),-u2(i),dpc)
enddo
end subroutine ztocd

!============================================================================
subroutine setmobius(xc0,xc1,xc2, aa,bb,cc,dd)!                   [setmobius]
!============================================================================
! Find the Mobius transformation complex coefficients, aa,bb,cc,dd,
! with aa*dd-bb*cc=1, for a standard (north-)polar stereographic transformation
! that takes cartesian point, xc0 to the north pole, xc1 to (lat=0,lon=0),
! xc2 to the south pole (=complex infinity).
!============================================================================
real(dp),dimension(3),intent(IN ):: xc0,xc1,xc2
complex(dpc),         intent(OUT):: aa,bb,cc,dd
!----------------------------------------------------------------------------
real(dp),parameter:: zero=0,one=1
logical           :: infz0,infz1,infz2
complex(dpc)      :: z0,z1,z2,z02,z10,z21
!============================================================================
call ctoz(xc0,z0,infz0)
call ctoz(xc1,z1,infz1)
call ctoz(xc2,z2,infz2)
z21=z2-z1
z02=z0-z2
z10=z1-z0

if(  (z0==z1.and.infz0.eqv.infz1).or.&
     (z1==z2.and.infz1.eqv.infz2).or.&
     (z2==z0.and.infz2.eqv.infz0))   &
     stop 'In setmobius; anchor points must be distinct'

if(infz2 .or. (.not.infz0 .and. abs(z0)<abs(z2)))then
! z0 is finite and smaller than z2:
   if(infz1)then
      aa=one/sqrt(z02)        ! <- z1 is infinite
   elseif(infz2)then
      aa=one/sqrt(z10)        ! <- z2 is infinite
   else
      aa=sqrt(-z21/(z02*z10)) ! <- all zs are finite
   endif
   bb=-z0*aa
   if(infz1)then
      cc=aa                   ! <- z1 is infinite
      dd=-z2*aa               !
   elseif(infz2)then
      cc=zero                 ! <- z2 is infinite
      dd=z10*aa               !
   else
      cc=-(z10/z21)*aa        ! <- all zs are finite
      dd= z2*(z10/z21)*aa     !
   endif
else
! z2 is finite and smaller than z0:
   if(infz0)then
      cc=one/sqrt(z21)        ! <- z0 is inifinite
   elseif(infz1)then
      cc=one/sqrt(z02)        ! <- z1 is infinite
   else
      cc=sqrt(-z10/(z02*z21)) ! <- all zs are finite
   endif
   dd=-z2*cc
   if(infz0)then
      aa=zero                 ! <- z0 is inifinite
      bb=-z21*cc              !
   elseif(infz1)then
      aa=cc                   ! <- z1 is infinite
      bb=-z0*cc               !
   else
      aa=(-z21/z10)*cc        ! <- all zs are finite
      bb=z0*(z21/z10)*cc      !
   endif
endif
end subroutine setmobius

!============================================================================
subroutine zsetmobius(z0,infz0, z1,infz1, z2,infz2,  aa,bb,cc,dd)
!                                                                 [setmobius]
!============================================================================
! Find the Mobius transformation complex coefficients, aa,bb,cc,dd,
! with aa*dd-bb*cc=1, 
! that takes polar stereographic  point, z0 to the north pole, 
! z1 to (lat=0,lon=0), z2 to the south pole (=complex infinity).
! Should any one of z0,z1,z2 be itself the "point at infinity" its 
! corresponding infz will be set "true" (and the z value itself not used).
! This routine is like setmobius, except the three fixed points defining
! the mapping are given in standard complex stereographic form, together
! with the logical codes "infzn" that are TRUE if that point is itself
! the projection pole (i.e., the South Pole for a north polar stereographic).
!============================================================================
complex(dp),          intent(IN ):: z0,z1,z2
logical,              intent(IN ):: infz0,infz1,infz2
complex(dpc),         intent(OUT):: aa,bb,cc,dd
!----------------------------------------------------------------------------
real(dp),parameter:: zero=0,one=1
complex(dpc)      :: z02,z10,z21
!============================================================================
z21=z2-z1
z02=z0-z2
z10=z1-z0

if(  (z0==z1.and.infz0.eqv.infz1).or.&
     (z1==z2.and.infz1.eqv.infz2).or.&
     (z2==z0.and.infz2.eqv.infz0))   &
     stop 'In setmobius; anchor points must be distinct'

if(infz2 .or. (.not.infz0 .and. abs(z0)<abs(z2)))then
! z0 is finite and smaller than z2:
   if(infz1)then
      aa=one/sqrt(z02)        ! <- z1 is infinite
   elseif(infz2)then
      aa=one/sqrt(z10)        ! <- z2 is infinite
   else
      aa=sqrt(-z21/(z02*z10)) ! <- all zs are finite
   endif
   bb=-z0*aa
   if(infz1)then
      cc=aa                   ! <- z1 is infinite
      dd=-z2*aa               !
   elseif(infz2)then
      cc=zero                 ! <- z2 is infinite
      dd=z10*aa               !
   else
      cc=-(z10/z21)*aa        ! <- all zs are finite
      dd= z2*(z10/z21)*aa     !
   endif
else
! z2 is finite and smaller than z0:
   if(infz0)then
      cc=one/sqrt(z21)        ! <- z0 is inifinite
   elseif(infz1)then
      cc=one/sqrt(z02)        ! <- z1 is infinite
   else
      cc=sqrt(-z10/(z02*z21)) ! <- all zs are finite
   endif
   dd=-z2*cc
   if(infz0)then
      aa=zero                 ! <- z0 is inifinite
      bb=-z21*cc              !
   elseif(infz1)then
      aa=cc                   ! <- z1 is infinite
      bb=-z0*cc               !
   else
      aa=(-z21/z10)*cc        ! <- all zs are finite
      bb=z0*(z21/z10)*cc      !
   endif
endif
end subroutine zsetmobius

!=============================================================================
subroutine zmobius(aa,bb,cc,dd, z,infz, w,infw)!                      [mobius]
!=============================================================================
! Perform a complex Mobius transformation from (z,infz) to (w,infw)
! where the transformation coefficients are the standard aa,bb,cc,dd.
! Infz is .TRUE. only when z is at complex infinity; likewise infw and w.
! For these infinite cases, it is important that numerical z==(0,0).
!=============================================================================
complex(dpc),intent(IN ):: aa,bb,cc,dd,z
logical,     intent(IN ):: infz
complex(dpc),intent(OUT):: w
logical,     intent(OUT):: infw
!-----------------------------------------------------------------------------
real(dp),parameter:: zero=0
complex(dpc)      :: top,bot
!=============================================================================
w=0
infw=.false.
if(infz)then
   top=aa
   bot=cc
else
   top=aa*z+bb
   bot=cc*z+dd
endif

if(abs(bot)==zero)then
   infw=.true.
else
   w=top/bot
endif
end subroutine zmobius

!=============================================================================
subroutine cmobius(aa,bb,cc,dd, vz,vw)!                               [mobius]
!=============================================================================
! Perform a complex Mobius transformation from cartesian vz to cartesian vw
! where the transformation coefficients are the standard aa,bb,cc,dd.
!=============================================================================
complex(dpc),         intent(IN ):: aa,bb,cc,dd
real(dp),dimension(3),intent(IN ):: vz
real(dp),dimension(3),intent(OUT):: vw
!-----------------------------------------------------------------------------
complex(dpc):: z,w
logical     :: infz,infw
!=============================================================================
call ctoz(vz, z,infz)
call zmobius(aa,bb,cc,dd, z,infz, w,infw)
call ztoc(w,infw, vw)
end subroutine cmobius

!============================================================================
subroutine zmobiusi(aa,bb,cc,dd, zz,infz, zw,infw)                ! [mobiusi]
!============================================================================
! Perform the inverse of the mobius transformation with coefficients,
! {aa,bb,cc,dd}.
!============================================================================
complex(dpc),intent(IN ):: aa,bb,cc,dd,zz
logical,     intent(IN ):: infz
complex(dpc),intent(OUT):: zw
logical,     intent(OUT):: infw
!----------------------------------------------------------------------------
real(dp),parameter:: one=1
complex(dpc)      :: aai,bbi,cci,ddi,d
!============================================================================
d=one/(aa*dd-bb*cc)
aai=dd*d
ddi=aa*d
bbi=-bb*d
cci=-cc*d
call zmobius(aai,bbi,cci,ddi, zz,infz, zw,infw)
end subroutine zmobiusi


end module peuc