LAPACK/ScaLAPACK Development

Hello,
I am new in using LAPACK subroutines. I need to use sgeqp3.f and sorgqr.f.
Before putting them in the fluid dynamics code I use, I made a trial, but I wasn't successful. Could you please help me in finding the mistake I am making? I think I didn't defined well some arguments like WORK and JPVT. How should I define them? The program runs with no mistakes but the TAU is never updated.
Down here there is the trial code
Thank you.
Chiara

Code: Select all

     program main



      integer m,n,LWORK,INFO
   integer,allocatable,dimension (:) :: JPVT
      real xb,yb
      real, allocatable, dimension (:) :: xbi,ybi,WORK,TAU
   real,allocatable,dimension (:,:) :: a   


      m=8
   n=5
   
   xb=2.
   yb=2.

      allocate(xbi(m))
   xbi(1)=1.
   xbi(2)=2.
   xbi(3)=3.
   xbi(4)=1.
   xbi(5)=3.
   xbi(6)=1.
   xbi(7)=2.
   xbi(8)=3.

      allocate(ybi(m))
   ybi(1)=2.
   ybi(2)=2.
   ybi(3)=2.
   ybi(4)=1.
   ybi(5)=1.
   ybi(6)=0.
   ybi(7)=0.
   ybi(8)=0.

      allocate(a(m,n))
   do i=1,m
   a(i,1)=xbi(i)-xb
   a(i,2)=ybi(i)-yb
   a(i,3)=(xbi(i)-xb)*(ybi(i)-yb)
   a(i,4)=0.5*(xbi(i)-xb)*(xbi(i)-xb)
      a(i,5)=0-5*(ybi(i)-yb)*(ybi(i)-yb)
   end do



   LWORK=-1
   allocate(WORK(1))
   allocate(TAU(n))
   allocate(JPVT(n))
   JPVT=0
      TAU=0.

   call SGEQP3(m,n,a,m,JPVT,TAU,WORK,LWORK,INFO)
!      LWORK=WORK(1)
!   deallocate(WORK)
!   allocate(WORK(LWORK))   
!   WORK=0
      write(*,*)JPVT(1:n)
   call SORGQR(m,n,n,a,m,TAU,WORK,LWORK,INFO)
     
      do i=1,m
   write(*,*)a(i,1:n)
   end do
   write(*,*)'WORK=',WORK(1)
   write(*,*)'INFO=',INFO
      write(*,*)

   do i=1,n
      write(*,*)TAU(i)
   end do


   end program

You might try reproducing the solution of a documented problem to make sure you are definiing things correctly:

http://www.nag.co.uk/numeric/fl/manual/ ... f08aff.pdf

Sorry, I think my post was lost somewhere in the system because I can not
find it anymore ...

Chiara, you need to call twice each LAPACK function: the first time you set
LWORK=-1 and this is a workspace query and you will get LWORK in
WORK(1), then you allocate WORK and then you call again the same LAPACK
function again and only now something interesting will really happen in
LAPACK.
So you need to call SGEQP3 twice, SORGQR twice.
Otherwise the code looks good.

BTW why do you use SGEQP3 and not SGEQRF?

Ciao,
thank you for your suggestions. I modified the code and try to use sgeqrf, but it still gives me problems since WORK(1) at the first passage returns 0. What is still wrong?
Thanks again.
Chiara

Code: Select all

     program main

      integer m,n,LWORK,INFO
      real, allocatable, dimension (:) :: WORK,TAU
   real,allocatable,dimension (:,:) :: a,b,c
   intrinsic dot_product
 
        m=4
   n=3
   
      allocate(a(m,n))
      allocate(b(n,m))
      allocate(c(m,m))

        a(1,1)=1.
   a(1,2)=1.
   a(1,3)=0.
   a(2,1)=-1.
   a(2,2)=0.
   a(2,3)=1.
   a(3,1)=0.
   a(3,2)=1.
   a(3,3)=1.
   a(4,1)=0.
   a(4,2)=0.
   a(4,3)=1.
   do i=1,m
          write(*,*)a(i,1:n)
   end do

   LWORK=-1
   allocate(WORK(1))
   allocate(TAU(n))
       TAU=0.

   call sgeqrf(m,n,a,m,TAU,WORK,LWORK,INFO)

       LWORK=WORK(1)
       deallocate(WORK)
       allocate(WORK(LWORK))   

       call sgeqrf(m,n,a,m,TAU,WORK,LWORK,INFO)

       deallocate(WORK)
 
        LWORK=-1
   allocate(WORK(1))
   call sorgqr(m,n,n,a,m,TAU,WORK,LWORK,INFO)

       LWORK=WORK(1)
       deallocate(WORK)
       allocate(WORK(LWORK))

       call sorgqr(m,n,n,a,m,TAU,WORK,LWORK,INFO)

   do i=1,m
         do j=1,n
        b(j,i)=a(i,j)
      end do
        end do
     
      c=0.
   call proma(m,m,n,a,b,c)
   write(*,*)
      write(*,*)'Q*Qt'
   do i=1,m
     write(*,*)c(i,1:m)
   end do

   end program

All you are doing in your LAPACK's portion of the code is correct and makes
sense. In input A is an m-by-n input matrix after the last sorgqr call, A
contains an orthonormal basis for the initial column set.
You can check that for example Q^T Q = I. In this equality, Q is indeed
the m-by-n output A of LAPACK. I is the n-by-n identity matrix.
Julien.