define(s_type,`real(skind)')
define(d_type,`real(dkind)')
define(c_type,`complex(skind)')
define(z_type,`complex(dkind)')
define(i_type,`integer')
define(TYPE,`$1_type')

module lasvd
implicit none
private
integer,parameter:: skind = kind(1.e0), dkind = kind(1.d0)

! sizeof is a less common extension, so provide sizes manually
integer,parameter:: i_size = 4
integer,parameter:: s_size = 4 
integer,parameter:: d_size = 8 
integer,parameter:: c_size = 2*s_size
integer,parameter:: z_size = 2*d_size

s_type,target:: s_empty
d_type,target:: d_empty
c_type,target:: c_empty
z_type,target:: z_empty

ifdef(`EXTERNALLOC',`integer,external:: loc')
ifdef(`UNSAFE',`! this is unsafe version',`define(SAFE)
! this is safe version')

define(MAKEPUBLIC,`
ifdef(``ISPUBLIC'$1`'',`! $1 already public',`
public $1
define(``ISPUBLIC'$1`'')
')
')

define(BIND,`interface $1
  module procedure s$2_95
  module procedure d$2_95
  module procedure c$2_95
  module procedure z$2_95
end interface
MAKEPUBLIC($1)')
define(BINDR,`interface $1
  module procedure s$2_95
  module procedure d$2_95
end interface
MAKEPUBLIC($1)')
define(BINDC,`interface $1
  module procedure c$2_95
  module procedure z$2_95
end interface
MAKEPUBLIC($1)')

BIND(GESVD,GESVD)
BIND(GESDD,GESDD)

contains

ifdef(`SAFE',`
define(VECEXT,`
subroutine $1vecext(vec,ppvec,wwvec)
$1_type,target:: vec(:)
$1_type,pointer:: ppvec,wwvec(:)
if (size(vec) > 1) then
  if ((loc(vec(2)) - loc(vec(1))) /= $1_size) then
    allocate(wwvec(size(vec)))
    ppvec => wwvec(1)
    return
  end if
end if
ppvec => vec(1)
end subroutine
')
VECEXT(s)
VECEXT(d)
VECEXT(c)
VECEXT(z)
VECEXT(i)

define(MATEXT,`
subroutine $1matext(mat,ppmat,ldmat,wwmat)
$1_type,target:: mat(:,:)
$1_type,pointer:: ppmat,wwmat(:,:)
integer,intent(out):: ldmat
if (size(mat,1) > 1) then
  if (loc(mat(2,1)) - loc(mat(1,1)) /= $1_size) goto 1
end if
ldmat = 1
if (size(mat,2) > 1) ldmat = (loc(mat(1,2)) - loc(mat(1,1))) / $1_size
if (ldmat < 0) goto 1
ppmat => mat(1,1)
return
1 continue
allocate(wwmat(size(mat,1),size(mat,2)))
ppmat => wwmat(1,1)
ldmat = size(mat,1)
end subroutine
')
MATEXT(s)
MATEXT(d)
MATEXT(c)
MATEXT(z)
')

ifdef(`SAFE',`
define(VECASS,`associated(pp$1,$1(1))')
define(VECARG,`
$3_type,dimension(:),intent($2),target:: $1
$3_type,pointer:: pp$1,ww$1(:)
pushdef(`COPY',`
call $3vecext($1,pp$1,ww$1)
ifelse($2,out,,`if (.not.VECASS($1)) ww$1 = $1')
popdef(`COPY')
COPY
')
pushdef(`CLEANUP',`
if (.not.VECASS($1)) then
  ifelse($2,in,,`$1 = ww$1')
  deallocate(pp$1)
end if
popdef(`CLEANUP')
CLEANUP
')
')
',`
define(VECARG,`
$3_type,dimension(:),intent($2),target:: $1
$3_type,pointer:: pp$1
pushdef(`COPY',`
pp$1 => $1(1)
popdef(`COPY')
COPY
')
')
')

ifdef(`SAFE',`
define(MATASS,`associated(pp$1,$1(1,1))')
define(MATARG,`
TYPE(XC),dimension(:,:),intent($2),target:: $1
TYPE(XC),pointer:: pp$1,ww$1(:,:)
integer:: ld$1
pushdef(`COPY',`
call XC`matext'($1,pp$1,ld$1,ww$1)
ifelse($2,out,,`if (.not.MATASS($1)) ww$1 = $1')
popdef(`COPY')
COPY
')
pushdef(`CLEANUP',`
if (.not.MATASS($1)) then
  ifelse($2,in,,`ifelse($3,,`$1 = ww$1',`if ($3) $1 = ww$1')')
  deallocate(pp$1)
end if
popdef(`CLEANUP')
CLEANUP
')
')
',`
define(MATARG,`
TYPE(XC),dimension(:,:),intent($2),target:: $1
TYPE(XC),pointer:: pp$1
integer:: ld$1
pushdef(`COPY',`
ld$1 = 1
pp$1 => $1(1,1)
if (size($1,2) > 1) ld$1 = (loc($1(1,2)) - loc($1(1,1)))/XC`_size'
popdef(`COPY')
COPY
')
')
')

ifdef(`SAFE',`
define(MATARGOO,`
TYPE(XC),dimension(:,:),intent(out),target,optional:: $1
TYPE(XC),pointer:: pp$1,ww$1(:,:)
integer:: ld$1
pushdef(`COPY',`
if (present($1)) then
  call XC`matext'($1,pp$1,ld$1,ww$1)
  if (.not.MATASS($1)) ww$1 = $1
else
  pp$1 => XC`_empty'
  ld$1 = 1
end if
popdef(`COPY')
COPY
')
pushdef(`CLEANUP',`
if (present($1)) then
  if (.not.MATASS($1)) then
    $1 = ww$1
    deallocate(pp$1)
  end if
end if
popdef(`CLEANUP')
CLEANUP
')
')
',`
define(MATARGOO,`
TYPE(XC),dimension(:,:),intent(out),target,optional:: $1
TYPE(XC),pointer:: pp$1
integer:: ld$1
pushdef(`COPY',`
if (present($1)) then
  ld$1 = 1
  pp$1 => $1(1,1)
  if (size($1,2) > 1) ld$1 = (loc($1(1,2)) - loc($1(1,1)))/XC`_size'
else
  pp$1 => XC`_empty'
  ld$1 = 1
end if
popdef(`COPY')
COPY
')
')
')

define(INFOARG,`
integer,intent(out),optional:: info
integer:: vinfo
pushdef(`CLEANUP',`
if (present(info)) info = vinfo
popdef(`CLEANUP')
CLEANUP
')
')


define(CHARG,`character,intent(in),optional:: $1
character:: v$1
pushdef(`INIT',`
v$1 = $2
if (present($1)) v$1 = $1
popdef(`INIT')
INIT
')
')
define(WORK,`
TYPE(XC),dimension(:),allocatable:: work
TYPE(XC):: lwork
')

define(SIZE,`
integer:: $1
pushdef(`INIT',`
$1 = size($2,$3)
if ($1 == 0) return
popdef(`INIT')
INIT
')
')

define(`INIT',`! end init')
define(`COPY',`! end copy')
define(`CLEANUP',`! end cleanup')

! singular value decompositions

define(GESVD,`
pushdef(`XC',$1)
subroutine $1GESVD_95(a,s,u,vt,ww,job,info)
MATARG(a,inout,`vjob == "N"')
VECARG(s,inout,$1)
$1_type,dimension(:),intent(out),optional:: ww
MATARGOO(u)
MATARGOO(vt)
CHARG(job,"N")
INFOARG()
SIZE(m,a,1)
SIZE(n,a,2)
WORK()
character:: jobu,jobvt
INIT
COPY
if (present(u)) then
  jobu = "S"
  if (size(u,2) == m) jobu = "A"
else
  jobu = "N"
  if (vjob == "U") jobu = "0"
end if
if (present(vt)) then
  jobvt = "S"
  if (size(vt,1) == n) jobvt = "A"
else
  jobvt = "N"
  if (vjob == "V") jobvt = "0"
end if
call $1GESVD(jobu,jobvt,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1GESVD(jobu,jobvt,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,work,size(work),vinfo)
if (vinfo == 0 .and. present(ww)) ww = work(2:min(m,n))
CLEANUP
popdef(`XC')
end subroutine
pushdef(`XC',$2)
subroutine $2GESVD_95(a,s,u,vt,ww,job,info)
MATARG(a,inout,`vjob == "N"')
VECARG(s,inout,$1)
$1_type,dimension(:),intent(out),optional:: ww
MATARGOO(u)
MATARGOO(vt)
CHARG(job,"N")
INFOARG()
SIZE(m,a,1)
SIZE(n,a,2)
WORK()
$1_type,dimension(:),allocatable:: rwork
character:: jobu,jobvt
INIT
COPY
if (present(u)) then
  jobu = "S"
  if (size(u,2) == m) jobu = "A"
else
  jobu = "N"
  if (vjob == "U") jobu = "0"
end if
if (present(vt)) then
  jobvt = "S"
  if (size(vt,1) == n) jobvt = "A"
else
  jobvt = "N"
  if (vjob == "V") jobvt = "0"
end if
allocate(rwork(5*min(m,n)))
call $2GESVD(jobu,jobvt,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,&
lwork,-1,rwork,vinfo)
allocate(work(int(lwork)))
call $2GESVD(jobu,jobvt,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,&
work,size(work),rwork,vinfo)
if (vinfo == 0 .and. present(ww)) ww = rwork(1:min(m,n)-1)
CLEANUP
popdef(`XC')
end subroutine
')
GESVD(s,c)
GESVD(d,z)

define(GESDD,`
pushdef(`XC',$1)
subroutine $1GESDD_95(a,s,u,vt,info)
MATARG(a,inout,`vjob == "N"')
VECARG(s,inout,$1)
MATARGOO(u)
MATARGOO(vt)
INFOARG()
SIZE(m,a,1)
SIZE(n,a,2)
WORK()
integer,dimension(:),allocatable:: iwork
character:: jobz
INIT
COPY
if (present(u).and.present(vt)) then
  jobz = "S"
  if (size(u,2) == m .and. size(vt,1) == n) jobz = "A"
else
  jobz = "0"
  if (.not.present(u).and..not.present(vt)) jobz = "N"
end if
allocate(iwork(8*min(m,n)))
call $1GESDD(jobz,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,lwork,-1,iwork,vinfo)
allocate(work(int(lwork)))
call $1GESDD(jobz,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,work,size(work),iwork,vinfo)
CLEANUP
popdef(`XC')
end subroutine
pushdef(`XC',$2)
subroutine $2GESDD_95(a,s,u,vt,job,info)
MATARG(a,inout,`vjob == "N"')
VECARG(s,inout,$1)
MATARGOO(u)
MATARGOO(vt)
CHARG(job,"N")
INFOARG()
SIZE(m,a,1)
SIZE(n,a,2)
WORK()
$1_type,dimension(:),allocatable:: rwork
integer,dimension(:),allocatable:: iwork
character:: jobz
INIT
COPY
if (present(u).and.present(vt)) then
  jobz = "S"
  if (size(u,2) == m .and. size(vt,1) == n) jobz = "A"
else
  jobz = "0"
  if (.not.present(u).and..not.present(vt)) jobz = "N"
end if
allocate(iwork(8*min(m,n)))
allocate(rwork(5*min(m,n)))
call $2GESDD(jobz,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,&
lwork,-1,rwork,iwork,vinfo)
allocate(work(int(lwork)))
call $2GESDD(jobz,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,&
work,size(work),rwork,iwork,vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
GESDD(s,c)
GESDD(d,z)

end module