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 lalsq
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(GEQRF,GEQRF)
BINDR(ORMQR,ORMQR)
BINDC(UNMQR,UNMQR)
BINDR(ORGQR,ORGQR)
BINDC(UNGQR,UNGQR)

BIND(GELQF,GELQF)
BINDR(ORMLQ,ORMLQ)
BINDC(UNMLQ,UNMLQ)
BINDR(ORGLQ,ORGLQ)
BINDC(UNGLQ,UNGLQ)

BIND(GEQP3,GEQP3)

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
')
')
')


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')

! orthogonal factorizations

! QR factorization
define(GEQRF,`
pushdef(`XC',$1)
subroutine $1GEQRF_95(a,tau,info)
MATARG(a,inout)
VECARG(tau,out,$1)
INFOARG()
WORK()
SIZE(m,a,1)
SIZE(n,a,2)
INIT
COPY
call $1GEQRF(m,n,ppa,lda,pptau,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1GEQRF(m,n,ppa,lda,pptau,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
GEQRF(s)
GEQRF(d)
GEQRF(c)
GEQRF(z)

define(ORMQR,`
pushdef(`XC',$1)
subroutine $1ORMQR_95(a,tau,c,side,trans,info)
MATARG(a,in)
VECARG(tau,in,$1)
MATARG(c,inout)
CHARG(side,"L")
CHARG(trans,"N")
INFOARG()
WORK()
INIT
COPY
call $1ORMQR(vside,vtrans,size(c,1),size(c,2),size(tau),&
ppa,lda,pptau,ppc,ldc,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1ORMQR(vside,vtrans,size(c,1),size(c,2),size(tau),&
ppa,lda,pptau,ppc,ldc,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
ORMQR(s)
ORMQR(d)

define(UNMQR,`
pushdef(`XC',$1)
subroutine $1UNMQR_95(a,tau,c,side,trans,info)
MATARG(a,in)
VECARG(tau,in,$1)
MATARG(c,inout)
CHARG(side,"L")
CHARG(trans,"N")
INFOARG()
WORK()
INIT
COPY
call $1UNMQR(vside,vtrans,size(c,1),size(c,2),size(tau),&
ppa,lda,pptau,ppc,ldc,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1UNMQR(vside,vtrans,size(c,1),size(c,2),size(tau),&
ppa,lda,pptau,ppc,ldc,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
UNMQR(c)
UNMQR(z)

define(ORGQR,`
pushdef(`XC',$1)
subroutine $1ORGQR_95(a,tau,info)
MATARG(a,inout)
VECARG(tau,in,$1)
INFOARG()
WORK()
INIT
COPY
call $1ORGQR(size(a,1),size(a,2),size(tau),&
ppa,lda,pptau,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1ORGQR(size(a,1),size(a,2),size(tau),&
ppa,lda,pptau,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
ORGQR(s)
ORGQR(d)

define(UNGQR,`
pushdef(`XC',$1)
subroutine $1UNGQR_95(a,tau,info)
MATARG(a,inout)
VECARG(tau,in,$1)
INFOARG()
WORK()
INIT
COPY
call $1UNGQR(size(a,1),size(a,2),size(tau),&
ppa,lda,pptau,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1UNGQR(size(a,1),size(a,2),size(tau),&
ppa,lda,pptau,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
UNGQR(c)
UNGQR(z)

! RQ factorization
define(GERQF,`
pushdef(`XC',$1)
subroutine $1GERQF_95(a,tau,info)
MATARG(a,inout)
VECARG(tau,out,$1)
INFOARG()
WORK()
SIZE(m,a,1)
SIZE(n,a,2)
INIT
COPY
call $1GERQF(m,n,ppa,lda,pptau,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1GERQF(m,n,ppa,lda,pptau,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
GERQF(s)
GERQF(d)
GERQF(c)
GERQF(z)

define(ORMRQ,`
pushdef(`XC',$1)
subroutine $1ORMRQ_95(a,tau,c,side,trans,info)
MATARG(a,in)
VECARG(tau,in,$1)
MATARG(c,inout)
CHARG(side,"L")
CHARG(trans,"N")
INFOARG()
WORK()
INIT
COPY
call $1ORMRQ(vside,vtrans,size(c,1),size(c,2),size(tau),&
ppa,lda,pptau,ppc,ldc,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1ORMRQ(vside,vtrans,size(c,1),size(c,2),size(tau),&
ppa,lda,pptau,ppc,ldc,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
ORMRQ(s)
ORMRQ(d)

define(UNMRQ,`
pushdef(`XC',$1)
subroutine $1UNMRQ_95(a,tau,c,side,trans,info)
MATARG(a,in)
VECARG(tau,in,$1)
MATARG(c,inout)
CHARG(side,"L")
CHARG(trans,"N")
INFOARG()
WORK()
INIT
COPY
call $1UNMRQ(vside,vtrans,size(c,1),size(c,2),size(tau),&
ppa,lda,pptau,ppc,ldc,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1UNMRQ(vside,vtrans,size(c,1),size(c,2),size(tau),&
ppa,lda,pptau,ppc,ldc,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
UNMRQ(c)
UNMRQ(z)

define(ORGRQ,`
pushdef(`XC',$1)
subroutine $1ORGRQ_95(a,tau,info)
MATARG(a,inout)
VECARG(tau,in,$1)
INFOARG()
WORK()
INIT
COPY
call $1ORGRQ(size(a,1),size(a,2),size(tau),&
ppa,lda,pptau,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1ORGRQ(size(a,1),size(a,2),size(tau),&
ppa,lda,pptau,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
ORGRQ(s)
ORGRQ(d)

define(UNGRQ,`
pushdef(`XC',$1)
subroutine $1UNGRQ_95(a,tau,info)
MATARG(a,inout)
VECARG(tau,in,$1)
INFOARG()
WORK()
INIT
COPY
call $1UNGRQ(size(a,1),size(a,2),size(tau),&
ppa,lda,pptau,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1UNGRQ(size(a,1),size(a,2),size(tau),&
ppa,lda,pptau,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
UNGRQ(c)
UNGRQ(z)

! LQ factorization
define(GELQF,`
pushdef(`XC',$1)
subroutine $1GELQF_95(a,tau,info)
MATARG(a,inout)
VECARG(tau,out,$1)
INFOARG()
WORK()
SIZE(m,a,1)
SIZE(n,a,2)
INIT
COPY
call $1GELQF(m,n,ppa,lda,pptau,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1GELQF(m,n,ppa,lda,pptau,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
GELQF(s)
GELQF(d)
GELQF(c)
GELQF(z)

define(ORMLQ,`
pushdef(`XC',$1)
subroutine $1ORMLQ_95(a,tau,c,side,trans,info)
MATARG(a,in)
VECARG(tau,in,$1)
MATARG(c,inout)
CHARG(side,"L")
CHARG(trans,"N")
INFOARG()
WORK()
INIT
COPY
call $1ORMLQ(vside,vtrans,size(c,1),size(c,2),size(tau),&
ppa,lda,pptau,ppc,ldc,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1ORMLQ(vside,vtrans,size(c,1),size(c,2),size(tau),&
ppa,lda,pptau,ppc,ldc,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
ORMLQ(s)
ORMLQ(d)

define(UNMLQ,`
pushdef(`XC',$1)
subroutine $1UNMLQ_95(a,tau,c,side,trans,info)
MATARG(a,in)
VECARG(tau,in,$1)
MATARG(c,inout)
CHARG(side,"L")
CHARG(trans,"N")
INFOARG()
WORK()
INIT
COPY
call $1UNMLQ(vside,vtrans,size(c,1),size(c,2),size(tau),&
ppa,lda,pptau,ppc,ldc,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1UNMLQ(vside,vtrans,size(c,1),size(c,2),size(tau),&
ppa,lda,pptau,ppc,ldc,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
UNMLQ(c)
UNMLQ(z)

define(ORGLQ,`
pushdef(`XC',$1)
subroutine $1ORGLQ_95(a,tau,info)
MATARG(a,inout)
VECARG(tau,in,$1)
INFOARG()
WORK()
INIT
COPY
call $1ORGLQ(size(a,1),size(a,2),size(tau),&
ppa,lda,pptau,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1ORGLQ(size(a,1),size(a,2),size(tau),&
ppa,lda,pptau,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
ORGLQ(s)
ORGLQ(d)

define(UNGLQ,`
pushdef(`XC',$1)
subroutine $1UNGLQ_95(a,tau,info)
MATARG(a,inout)
VECARG(tau,in,$1)
INFOARG()
WORK()
INIT
COPY
call $1UNGLQ(size(a,1),size(a,2),size(tau),&
ppa,lda,pptau,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1UNGLQ(size(a,1),size(a,2),size(tau),&
ppa,lda,pptau,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
UNGLQ(c)
UNGLQ(z)

! QR factorization with column pivoting
define(GEQP3,`
pushdef(`XC',$1)
subroutine $1GEQP3_95(a,jpvt,tau,info)
MATARG(a,inout)
VECARG(jpvt,out,i)
VECARG(tau,out,$1)
INFOARG()
WORK()
INIT
COPY
call $1GEQP3(size(a,1),size(a,2),ppa,lda,&
ppjpvt,pptau,lwork,-1,vinfo)
allocate(work(int(lwork)))
call $1GEQP3(size(a,1),size(a,2),ppa,lda,&
ppjpvt,pptau,work,size(work),vinfo)
CLEANUP
popdef(`XC')
end subroutine
')
GEQP3(s)
GEQP3(d)
GEQP3(c)
GEQP3(z)

end module