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

real(skind),target:: s_empty
real(dkind),target:: d_empty
complex(skind),target:: c_empty
complex(dkind),target:: z_empty



! this is safe version







interface GESVD
  module procedure sGESVD_95
  module procedure dGESVD_95
  module procedure cGESVD_95
  module procedure zGESVD_95
end interface


public GESVD



interface GESDD
  module procedure sGESDD_95
  module procedure dGESDD_95
  module procedure cGESDD_95
  module procedure zGESDD_95
end interface


public GESDD




contains




subroutine svecext(vec,ppvec,wwvec)
real(skind),target:: vec(:)
real(skind),pointer:: ppvec,wwvec(:)
if (size(vec) > 1) then
  if ((loc(vec(2)) - loc(vec(1))) /= s_size) then
    allocate(wwvec(size(vec)))
    ppvec => wwvec(1)
    return
  end if
end if
ppvec => vec(1)
end subroutine


subroutine dvecext(vec,ppvec,wwvec)
real(dkind),target:: vec(:)
real(dkind),pointer:: ppvec,wwvec(:)
if (size(vec) > 1) then
  if ((loc(vec(2)) - loc(vec(1))) /= d_size) then
    allocate(wwvec(size(vec)))
    ppvec => wwvec(1)
    return
  end if
end if
ppvec => vec(1)
end subroutine


subroutine cvecext(vec,ppvec,wwvec)
complex(skind),target:: vec(:)
complex(skind),pointer:: ppvec,wwvec(:)
if (size(vec) > 1) then
  if ((loc(vec(2)) - loc(vec(1))) /= c_size) then
    allocate(wwvec(size(vec)))
    ppvec => wwvec(1)
    return
  end if
end if
ppvec => vec(1)
end subroutine


subroutine zvecext(vec,ppvec,wwvec)
complex(dkind),target:: vec(:)
complex(dkind),pointer:: ppvec,wwvec(:)
if (size(vec) > 1) then
  if ((loc(vec(2)) - loc(vec(1))) /= z_size) then
    allocate(wwvec(size(vec)))
    ppvec => wwvec(1)
    return
  end if
end if
ppvec => vec(1)
end subroutine


subroutine ivecext(vec,ppvec,wwvec)
integer,target:: vec(:)
integer,pointer:: ppvec,wwvec(:)
if (size(vec) > 1) then
  if ((loc(vec(2)) - loc(vec(1))) /= i_size) then
    allocate(wwvec(size(vec)))
    ppvec => wwvec(1)
    return
  end if
end if
ppvec => vec(1)
end subroutine




subroutine smatext(mat,ppmat,ldmat,wwmat)
real(skind),target:: mat(:,:)
real(skind),pointer:: ppmat,wwmat(:,:)
integer,intent(out):: ldmat
if (size(mat,1) > 1) then
  if (loc(mat(2,1)) - loc(mat(1,1)) /= s_size) goto 1
end if
ldmat = 1
if (size(mat,2) > 1) ldmat = (loc(mat(1,2)) - loc(mat(1,1))) / s_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


subroutine dmatext(mat,ppmat,ldmat,wwmat)
real(dkind),target:: mat(:,:)
real(dkind),pointer:: ppmat,wwmat(:,:)
integer,intent(out):: ldmat
if (size(mat,1) > 1) then
  if (loc(mat(2,1)) - loc(mat(1,1)) /= d_size) goto 1
end if
ldmat = 1
if (size(mat,2) > 1) ldmat = (loc(mat(1,2)) - loc(mat(1,1))) / d_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


subroutine cmatext(mat,ppmat,ldmat,wwmat)
complex(skind),target:: mat(:,:)
complex(skind),pointer:: ppmat,wwmat(:,:)
integer,intent(out):: ldmat
if (size(mat,1) > 1) then
  if (loc(mat(2,1)) - loc(mat(1,1)) /= c_size) goto 1
end if
ldmat = 1
if (size(mat,2) > 1) ldmat = (loc(mat(1,2)) - loc(mat(1,1))) / c_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


subroutine zmatext(mat,ppmat,ldmat,wwmat)
complex(dkind),target:: mat(:,:)
complex(dkind),pointer:: ppmat,wwmat(:,:)
integer,intent(out):: ldmat
if (size(mat,1) > 1) then
  if (loc(mat(2,1)) - loc(mat(1,1)) /= z_size) goto 1
end if
ldmat = 1
if (size(mat,2) > 1) ldmat = (loc(mat(1,2)) - loc(mat(1,1))) / z_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





























! singular value decompositions




subroutine sGESVD_95(a,s,u,vt,ww,job,info)

real(skind),dimension(:,:),intent(inout),target:: a
real(skind),pointer:: ppa,wwa(:,:)
integer:: lda




real(skind),dimension(:),intent(inout),target:: s
real(skind),pointer:: pps,wws(:)



real(skind),dimension(:),intent(out),optional:: ww

real(skind),dimension(:,:),intent(out),target,optional:: u
real(skind),pointer:: ppu,wwu(:,:)
integer:: ldu




real(skind),dimension(:,:),intent(out),target,optional:: vt
real(skind),pointer:: ppvt,wwvt(:,:)
integer:: ldvt



character,intent(in),optional:: job
character:: vjob



integer,intent(out),optional:: info
integer:: vinfo



integer:: m



integer:: n



real(skind),dimension(:),allocatable:: work
real(skind):: lwork

character:: jobu,jobvt

n = size(a,2)
if (n == 0) return


m = size(a,1)
if (m == 0) return


vjob = "N"
if (present(job)) vjob = job

! end init




if (present(vt)) then
  call smatext(vt,ppvt,ldvt,wwvt)
  if (.not.associated(ppvt,vt(1,1))) wwvt = vt
else
  ppvt => s_empty
  ldvt = 1
end if


if (present(u)) then
  call smatext(u,ppu,ldu,wwu)
  if (.not.associated(ppu,u(1,1))) wwu = u
else
  ppu => s_empty
  ldu = 1
end if


call svecext(s,pps,wws)
if (.not.associated(pps,s(1))) wws = s


call smatext(a,ppa,lda,wwa)
if (.not.associated(ppa,a(1,1))) wwa = a

! end 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 sGESVD(jobu,jobvt,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,lwork,-1,vinfo)
allocate(work(int(lwork)))
call sGESVD(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))

if (present(info)) info = vinfo


if (present(vt)) then
  if (.not.associated(ppvt,vt(1,1))) then
    vt = wwvt
    deallocate(ppvt)
  end if
end if


if (present(u)) then
  if (.not.associated(ppu,u(1,1))) then
    u = wwu
    deallocate(ppu)
  end if
end if


if (.not.associated(pps,s(1))) then
  s = wws
  deallocate(pps)
end if


if (.not.associated(ppa,a(1,1))) then
  if (vjob == "N") a = wwa
  deallocate(ppa)
end if

! end cleanup






end subroutine

subroutine cGESVD_95(a,s,u,vt,ww,job,info)

complex(skind),dimension(:,:),intent(inout),target:: a
complex(skind),pointer:: ppa,wwa(:,:)
integer:: lda




real(skind),dimension(:),intent(inout),target:: s
real(skind),pointer:: pps,wws(:)



real(skind),dimension(:),intent(out),optional:: ww

complex(skind),dimension(:,:),intent(out),target,optional:: u
complex(skind),pointer:: ppu,wwu(:,:)
integer:: ldu




complex(skind),dimension(:,:),intent(out),target,optional:: vt
complex(skind),pointer:: ppvt,wwvt(:,:)
integer:: ldvt



character,intent(in),optional:: job
character:: vjob



integer,intent(out),optional:: info
integer:: vinfo



integer:: m



integer:: n



complex(skind),dimension(:),allocatable:: work
complex(skind):: lwork

real(skind),dimension(:),allocatable:: rwork
character:: jobu,jobvt

n = size(a,2)
if (n == 0) return


m = size(a,1)
if (m == 0) return


vjob = "N"
if (present(job)) vjob = job

! end init




if (present(vt)) then
  call cmatext(vt,ppvt,ldvt,wwvt)
  if (.not.associated(ppvt,vt(1,1))) wwvt = vt
else
  ppvt => c_empty
  ldvt = 1
end if


if (present(u)) then
  call cmatext(u,ppu,ldu,wwu)
  if (.not.associated(ppu,u(1,1))) wwu = u
else
  ppu => c_empty
  ldu = 1
end if


call svecext(s,pps,wws)
if (.not.associated(pps,s(1))) wws = s


call cmatext(a,ppa,lda,wwa)
if (.not.associated(ppa,a(1,1))) wwa = a

! end 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 cGESVD(jobu,jobvt,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,&
lwork,-1,rwork,vinfo)
allocate(work(int(lwork)))
call cGESVD(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)

if (present(info)) info = vinfo


if (present(vt)) then
  if (.not.associated(ppvt,vt(1,1))) then
    vt = wwvt
    deallocate(ppvt)
  end if
end if


if (present(u)) then
  if (.not.associated(ppu,u(1,1))) then
    u = wwu
    deallocate(ppu)
  end if
end if


if (.not.associated(pps,s(1))) then
  s = wws
  deallocate(pps)
end if


if (.not.associated(ppa,a(1,1))) then
  if (vjob == "N") a = wwa
  deallocate(ppa)
end if

! end cleanup






end subroutine



subroutine dGESVD_95(a,s,u,vt,ww,job,info)

real(dkind),dimension(:,:),intent(inout),target:: a
real(dkind),pointer:: ppa,wwa(:,:)
integer:: lda




real(dkind),dimension(:),intent(inout),target:: s
real(dkind),pointer:: pps,wws(:)



real(dkind),dimension(:),intent(out),optional:: ww

real(dkind),dimension(:,:),intent(out),target,optional:: u
real(dkind),pointer:: ppu,wwu(:,:)
integer:: ldu




real(dkind),dimension(:,:),intent(out),target,optional:: vt
real(dkind),pointer:: ppvt,wwvt(:,:)
integer:: ldvt



character,intent(in),optional:: job
character:: vjob



integer,intent(out),optional:: info
integer:: vinfo



integer:: m



integer:: n



real(dkind),dimension(:),allocatable:: work
real(dkind):: lwork

character:: jobu,jobvt

n = size(a,2)
if (n == 0) return


m = size(a,1)
if (m == 0) return


vjob = "N"
if (present(job)) vjob = job

! end init




if (present(vt)) then
  call dmatext(vt,ppvt,ldvt,wwvt)
  if (.not.associated(ppvt,vt(1,1))) wwvt = vt
else
  ppvt => d_empty
  ldvt = 1
end if


if (present(u)) then
  call dmatext(u,ppu,ldu,wwu)
  if (.not.associated(ppu,u(1,1))) wwu = u
else
  ppu => d_empty
  ldu = 1
end if


call dvecext(s,pps,wws)
if (.not.associated(pps,s(1))) wws = s


call dmatext(a,ppa,lda,wwa)
if (.not.associated(ppa,a(1,1))) wwa = a

! end 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 dGESVD(jobu,jobvt,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,lwork,-1,vinfo)
allocate(work(int(lwork)))
call dGESVD(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))

if (present(info)) info = vinfo


if (present(vt)) then
  if (.not.associated(ppvt,vt(1,1))) then
    vt = wwvt
    deallocate(ppvt)
  end if
end if


if (present(u)) then
  if (.not.associated(ppu,u(1,1))) then
    u = wwu
    deallocate(ppu)
  end if
end if


if (.not.associated(pps,s(1))) then
  s = wws
  deallocate(pps)
end if


if (.not.associated(ppa,a(1,1))) then
  if (vjob == "N") a = wwa
  deallocate(ppa)
end if

! end cleanup






end subroutine

subroutine zGESVD_95(a,s,u,vt,ww,job,info)

complex(dkind),dimension(:,:),intent(inout),target:: a
complex(dkind),pointer:: ppa,wwa(:,:)
integer:: lda




real(dkind),dimension(:),intent(inout),target:: s
real(dkind),pointer:: pps,wws(:)



real(dkind),dimension(:),intent(out),optional:: ww

complex(dkind),dimension(:,:),intent(out),target,optional:: u
complex(dkind),pointer:: ppu,wwu(:,:)
integer:: ldu




complex(dkind),dimension(:,:),intent(out),target,optional:: vt
complex(dkind),pointer:: ppvt,wwvt(:,:)
integer:: ldvt



character,intent(in),optional:: job
character:: vjob



integer,intent(out),optional:: info
integer:: vinfo



integer:: m



integer:: n



complex(dkind),dimension(:),allocatable:: work
complex(dkind):: lwork

real(dkind),dimension(:),allocatable:: rwork
character:: jobu,jobvt

n = size(a,2)
if (n == 0) return


m = size(a,1)
if (m == 0) return


vjob = "N"
if (present(job)) vjob = job

! end init




if (present(vt)) then
  call zmatext(vt,ppvt,ldvt,wwvt)
  if (.not.associated(ppvt,vt(1,1))) wwvt = vt
else
  ppvt => z_empty
  ldvt = 1
end if


if (present(u)) then
  call zmatext(u,ppu,ldu,wwu)
  if (.not.associated(ppu,u(1,1))) wwu = u
else
  ppu => z_empty
  ldu = 1
end if


call dvecext(s,pps,wws)
if (.not.associated(pps,s(1))) wws = s


call zmatext(a,ppa,lda,wwa)
if (.not.associated(ppa,a(1,1))) wwa = a

! end 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 zGESVD(jobu,jobvt,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,&
lwork,-1,rwork,vinfo)
allocate(work(int(lwork)))
call zGESVD(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)

if (present(info)) info = vinfo


if (present(vt)) then
  if (.not.associated(ppvt,vt(1,1))) then
    vt = wwvt
    deallocate(ppvt)
  end if
end if


if (present(u)) then
  if (.not.associated(ppu,u(1,1))) then
    u = wwu
    deallocate(ppu)
  end if
end if


if (.not.associated(pps,s(1))) then
  s = wws
  deallocate(pps)
end if


if (.not.associated(ppa,a(1,1))) then
  if (vjob == "N") a = wwa
  deallocate(ppa)
end if

! end cleanup






end subroutine





subroutine sGESDD_95(a,s,u,vt,info)

real(skind),dimension(:,:),intent(inout),target:: a
real(skind),pointer:: ppa,wwa(:,:)
integer:: lda




real(skind),dimension(:),intent(inout),target:: s
real(skind),pointer:: pps,wws(:)




real(skind),dimension(:,:),intent(out),target,optional:: u
real(skind),pointer:: ppu,wwu(:,:)
integer:: ldu




real(skind),dimension(:,:),intent(out),target,optional:: vt
real(skind),pointer:: ppvt,wwvt(:,:)
integer:: ldvt




integer,intent(out),optional:: info
integer:: vinfo



integer:: m



integer:: n



real(skind),dimension(:),allocatable:: work
real(skind):: lwork

integer,dimension(:),allocatable:: iwork
character:: jobz

n = size(a,2)
if (n == 0) return


m = size(a,1)
if (m == 0) return

! end init



if (present(vt)) then
  call smatext(vt,ppvt,ldvt,wwvt)
  if (.not.associated(ppvt,vt(1,1))) wwvt = vt
else
  ppvt => s_empty
  ldvt = 1
end if


if (present(u)) then
  call smatext(u,ppu,ldu,wwu)
  if (.not.associated(ppu,u(1,1))) wwu = u
else
  ppu => s_empty
  ldu = 1
end if


call svecext(s,pps,wws)
if (.not.associated(pps,s(1))) wws = s


call smatext(a,ppa,lda,wwa)
if (.not.associated(ppa,a(1,1))) wwa = a

! end 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 sGESDD(jobz,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,lwork,-1,iwork,vinfo)
allocate(work(int(lwork)))
call sGESDD(jobz,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,work,size(work),iwork,vinfo)

if (present(info)) info = vinfo


if (present(vt)) then
  if (.not.associated(ppvt,vt(1,1))) then
    vt = wwvt
    deallocate(ppvt)
  end if
end if


if (present(u)) then
  if (.not.associated(ppu,u(1,1))) then
    u = wwu
    deallocate(ppu)
  end if
end if


if (.not.associated(pps,s(1))) then
  s = wws
  deallocate(pps)
end if


if (.not.associated(ppa,a(1,1))) then
  if (vjob == "N") a = wwa
  deallocate(ppa)
end if

! end cleanup






end subroutine

subroutine cGESDD_95(a,s,u,vt,job,info)

complex(skind),dimension(:,:),intent(inout),target:: a
complex(skind),pointer:: ppa,wwa(:,:)
integer:: lda




real(skind),dimension(:),intent(inout),target:: s
real(skind),pointer:: pps,wws(:)




complex(skind),dimension(:,:),intent(out),target,optional:: u
complex(skind),pointer:: ppu,wwu(:,:)
integer:: ldu




complex(skind),dimension(:,:),intent(out),target,optional:: vt
complex(skind),pointer:: ppvt,wwvt(:,:)
integer:: ldvt



character,intent(in),optional:: job
character:: vjob



integer,intent(out),optional:: info
integer:: vinfo



integer:: m



integer:: n



complex(skind),dimension(:),allocatable:: work
complex(skind):: lwork

real(skind),dimension(:),allocatable:: rwork
integer,dimension(:),allocatable:: iwork
character:: jobz

n = size(a,2)
if (n == 0) return


m = size(a,1)
if (m == 0) return


vjob = "N"
if (present(job)) vjob = job

! end init




if (present(vt)) then
  call cmatext(vt,ppvt,ldvt,wwvt)
  if (.not.associated(ppvt,vt(1,1))) wwvt = vt
else
  ppvt => c_empty
  ldvt = 1
end if


if (present(u)) then
  call cmatext(u,ppu,ldu,wwu)
  if (.not.associated(ppu,u(1,1))) wwu = u
else
  ppu => c_empty
  ldu = 1
end if


call svecext(s,pps,wws)
if (.not.associated(pps,s(1))) wws = s


call cmatext(a,ppa,lda,wwa)
if (.not.associated(ppa,a(1,1))) wwa = a

! end 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 cGESDD(jobz,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,&
lwork,-1,rwork,iwork,vinfo)
allocate(work(int(lwork)))
call cGESDD(jobz,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,&
work,size(work),rwork,iwork,vinfo)

if (present(info)) info = vinfo


if (present(vt)) then
  if (.not.associated(ppvt,vt(1,1))) then
    vt = wwvt
    deallocate(ppvt)
  end if
end if


if (present(u)) then
  if (.not.associated(ppu,u(1,1))) then
    u = wwu
    deallocate(ppu)
  end if
end if


if (.not.associated(pps,s(1))) then
  s = wws
  deallocate(pps)
end if


if (.not.associated(ppa,a(1,1))) then
  if (vjob == "N") a = wwa
  deallocate(ppa)
end if

! end cleanup






end subroutine



subroutine dGESDD_95(a,s,u,vt,info)

real(dkind),dimension(:,:),intent(inout),target:: a
real(dkind),pointer:: ppa,wwa(:,:)
integer:: lda




real(dkind),dimension(:),intent(inout),target:: s
real(dkind),pointer:: pps,wws(:)




real(dkind),dimension(:,:),intent(out),target,optional:: u
real(dkind),pointer:: ppu,wwu(:,:)
integer:: ldu




real(dkind),dimension(:,:),intent(out),target,optional:: vt
real(dkind),pointer:: ppvt,wwvt(:,:)
integer:: ldvt




integer,intent(out),optional:: info
integer:: vinfo



integer:: m



integer:: n



real(dkind),dimension(:),allocatable:: work
real(dkind):: lwork

integer,dimension(:),allocatable:: iwork
character:: jobz

n = size(a,2)
if (n == 0) return


m = size(a,1)
if (m == 0) return

! end init



if (present(vt)) then
  call dmatext(vt,ppvt,ldvt,wwvt)
  if (.not.associated(ppvt,vt(1,1))) wwvt = vt
else
  ppvt => d_empty
  ldvt = 1
end if


if (present(u)) then
  call dmatext(u,ppu,ldu,wwu)
  if (.not.associated(ppu,u(1,1))) wwu = u
else
  ppu => d_empty
  ldu = 1
end if


call dvecext(s,pps,wws)
if (.not.associated(pps,s(1))) wws = s


call dmatext(a,ppa,lda,wwa)
if (.not.associated(ppa,a(1,1))) wwa = a

! end 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 dGESDD(jobz,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,lwork,-1,iwork,vinfo)
allocate(work(int(lwork)))
call dGESDD(jobz,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,work,size(work),iwork,vinfo)

if (present(info)) info = vinfo


if (present(vt)) then
  if (.not.associated(ppvt,vt(1,1))) then
    vt = wwvt
    deallocate(ppvt)
  end if
end if


if (present(u)) then
  if (.not.associated(ppu,u(1,1))) then
    u = wwu
    deallocate(ppu)
  end if
end if


if (.not.associated(pps,s(1))) then
  s = wws
  deallocate(pps)
end if


if (.not.associated(ppa,a(1,1))) then
  if (vjob == "N") a = wwa
  deallocate(ppa)
end if

! end cleanup






end subroutine

subroutine zGESDD_95(a,s,u,vt,job,info)

complex(dkind),dimension(:,:),intent(inout),target:: a
complex(dkind),pointer:: ppa,wwa(:,:)
integer:: lda




real(dkind),dimension(:),intent(inout),target:: s
real(dkind),pointer:: pps,wws(:)




complex(dkind),dimension(:,:),intent(out),target,optional:: u
complex(dkind),pointer:: ppu,wwu(:,:)
integer:: ldu




complex(dkind),dimension(:,:),intent(out),target,optional:: vt
complex(dkind),pointer:: ppvt,wwvt(:,:)
integer:: ldvt



character,intent(in),optional:: job
character:: vjob



integer,intent(out),optional:: info
integer:: vinfo



integer:: m



integer:: n



complex(dkind),dimension(:),allocatable:: work
complex(dkind):: lwork

real(dkind),dimension(:),allocatable:: rwork
integer,dimension(:),allocatable:: iwork
character:: jobz

n = size(a,2)
if (n == 0) return


m = size(a,1)
if (m == 0) return


vjob = "N"
if (present(job)) vjob = job

! end init




if (present(vt)) then
  call zmatext(vt,ppvt,ldvt,wwvt)
  if (.not.associated(ppvt,vt(1,1))) wwvt = vt
else
  ppvt => z_empty
  ldvt = 1
end if


if (present(u)) then
  call zmatext(u,ppu,ldu,wwu)
  if (.not.associated(ppu,u(1,1))) wwu = u
else
  ppu => z_empty
  ldu = 1
end if


call dvecext(s,pps,wws)
if (.not.associated(pps,s(1))) wws = s


call zmatext(a,ppa,lda,wwa)
if (.not.associated(ppa,a(1,1))) wwa = a

! end 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 zGESDD(jobz,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,&
lwork,-1,rwork,iwork,vinfo)
allocate(work(int(lwork)))
call zGESDD(jobz,m,n,ppa,lda,pps,ppu,ldu,ppvt,ldvt,&
work,size(work),rwork,iwork,vinfo)

if (present(info)) info = vinfo


if (present(vt)) then
  if (.not.associated(ppvt,vt(1,1))) then
    vt = wwvt
    deallocate(ppvt)
  end if
end if


if (present(u)) then
  if (.not.associated(ppu,u(1,1))) then
    u = wwu
    deallocate(ppu)
  end if
end if


if (.not.associated(pps,s(1))) then
  s = wws
  deallocate(pps)
end if


if (.not.associated(ppa,a(1,1))) then
  if (vjob == "N") a = wwa
  deallocate(ppa)
end if

! end cleanup






end subroutine


end module