Problem with the rank of the matrix

[happykaren100]

Hey guys,
I am writing a program relating the traffic problem. In some parts, I need to calculate the rank of the given matrix (not a square matrix). May I ask any function can be called in fortran to do so? Or what else can I write?

Thanks!

 

[mikrom]

You can use LAPACK - see this thread:

http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=818

and here

http://www.cs.ucdavis.edu/~bai/publications/baidemmeletal06.pdf

on page 75-16 you have an example how to use LAPACK subroutine SGESVD for computing rank of an 8-by-5 matrix.

 

[mikrom]

Other example - for DGESVD is here:

http://software.intel.com/sites/pro...mkl_sa/11/mkl_lapack_examples/dgesvd_ex.f.htm

 

[happykaren100]

Thanks for your reply!
May I ask how can I call the DGESVD function / use LAPACK ?
It shows "unresolved external symbol _DGESVD@64" when the program ran.

I am writing the program by the Microsoft Developer Studio. In fact I cannot find this function on it.....

 

[mikrom]

Hi happykaren100,
Look for exaple at this thread:

http://www.tek-tips.com/viewthread.cfm?qid=1678628

First you have to download the LAPACK-DLL and then link it with your program. In the mentioned thread it's only shown with gfortran compiler. With MS Studio you probably use the Intel Compiler or any othe commercial one - I have no experience with it. But you could try it first with gfortan (whih is free) and then when it works try it with your compiler.

 

[mikrom]

Hi Happy Karen,

Here is the working example:

[COLOR=#a020f0]Program[/color] RankOfMatrix
  [COLOR=#0000ff]! Computing singular values of a matrix using LAPACK[/color]
  [COLOR=#2e8b57][b]Implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]

  [COLOR=#0000ff]! declarations[/color]
  [COLOR=#2e8b57][b]integer[/b][/color] :: i, m, n, lda, lds, lwork, lwmax, info, rank
  [COLOR=#0000ff]! lwmax = max(3*min(m. n) + max(m,n), 5*min(m,n))[/color]
  [COLOR=#0000ff]!       = max(3*5 + 8, 25) = 25[/color]
  [COLOR=#2e8b57][b]parameter[/b][/color](m [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]8[/color], n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]5[/color], lwmax[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]25[/color])
  [COLOR=#0000ff]! lda = max(1, m)[/color]
  [COLOR=#0000ff]! lds = max(1, m, n)[/color]
  [COLOR=#2e8b57][b]parameter[/b][/color](lda [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]8[/color], lds [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]8[/color])
  [COLOR=#2e8b57][b]double precision[/b][/color] :: A(lda,n), S(lds), U(m,m), VT(n), work(lwmax), eps
 
  [COLOR=#0000ff]! matrix A[/color]
  A([COLOR=#ff00ff]1[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]22[/color], [COLOR=#ff00ff]10[/color],  [COLOR=#ff00ff]2[/color],  [COLOR=#ff00ff]3[/color],  [COLOR=#ff00ff]7[/color][COLOR=#804040][b]/[/b][/color])
  A([COLOR=#ff00ff]2[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]14[/color],  [COLOR=#ff00ff]7[/color], [COLOR=#ff00ff]10[/color],  [COLOR=#ff00ff]0[/color],  [COLOR=#ff00ff]8[/color][COLOR=#804040][b]/[/b][/color])  
  A([COLOR=#ff00ff]3[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/-[/b][/color][COLOR=#ff00ff]1[/color], [COLOR=#ff00ff]13[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]11[/color],  [COLOR=#ff00ff]3[/color][COLOR=#804040][b]/[/b][/color])   
  A([COLOR=#ff00ff]4[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/-[/b][/color][COLOR=#ff00ff]3[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]2[/color], [COLOR=#ff00ff]13[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]2[/color],  [COLOR=#ff00ff]4[/color][COLOR=#804040][b]/[/b][/color])    
  A([COLOR=#ff00ff]5[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]9[/color],  [COLOR=#ff00ff]8[/color],  [COLOR=#ff00ff]1[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]2[/color],  [COLOR=#ff00ff]4[/color][COLOR=#804040][b]/[/b][/color])
  A([COLOR=#ff00ff]6[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]9[/color],  [COLOR=#ff00ff]1[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]7[/color],  [COLOR=#ff00ff]5[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color][COLOR=#804040][b]/[/b][/color])
  A([COLOR=#ff00ff]7[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]2[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]6[/color],  [COLOR=#ff00ff]6[/color],  [COLOR=#ff00ff]5[/color],  [COLOR=#ff00ff]1[/color][COLOR=#804040][b]/[/b][/color])
  A([COLOR=#ff00ff]8[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]4[/color],  [COLOR=#ff00ff]5[/color],  [COLOR=#ff00ff]0[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]2[/color],  [COLOR=#ff00ff]2[/color][COLOR=#804040][b]/[/b][/color])


  [COLOR=#0000ff]! workspace query: calculates the optimal size of the WORK array[/color]
  lwork [COLOR=#804040][b]=[/b][/color] [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]
  [COLOR=#008080]CALL[/color] DGESVD([COLOR=#ff00ff]'N'[/color], [COLOR=#ff00ff]'N'[/color], m, n, A, lda, S, U, m, VT, n, WORK, lwork, info)
  
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'lwmax   = '[/color], lwmax
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'work(1) = '[/color], work([COLOR=#ff00ff]1[/color])
  [COLOR=#0000ff]! compute workspace size[/color]
  lwork [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]min[/color](lwmax, [COLOR=#008080]int[/color](work([COLOR=#ff00ff]1[/color])))
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'lwork = '[/color], lwork

  [COLOR=#0000ff]! solve[/color]
  [COLOR=#008080]CALL[/color] DGESVD([COLOR=#ff00ff]'N'[/color], [COLOR=#ff00ff]'N'[/color], m, n, A, lda, S, U, m, VT, n, WORK, lwork, info)

  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'info = '[/color], info
  [COLOR=#804040][b]if[/b][/color] (info [COLOR=#804040][b].eq.[/b][/color] [COLOR=#ff00ff]0[/color]) [COLOR=#804040][b]then[/b][/color]
    [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Singular values were succesfully computed:'[/color]
    [COLOR=#0000ff]! print the solution x[/color]
    [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], lds
      [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]9[/color]) i, S(i)
    [COLOR=#804040][b]end do[/b][/color]  
  [COLOR=#804040][b]else[/b][/color]
    [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'* Error computing singula values!'[/color]
  [COLOR=#804040][b]end if[/b][/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])

  [COLOR=#0000ff]! The rank of the matrix is the number of singular values that are not zero[/color]
  eps [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]5e-10[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'epsilon = '[/color], eps
  rank [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
  [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], lds
    [COLOR=#804040][b]if[/b][/color] (S(i) [COLOR=#804040][b].gt.[/b][/color] eps) [COLOR=#804040][b]then[/b][/color]
      rank [COLOR=#804040][b]=[/b][/color] rank [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
    [COLOR=#804040][b]end if[/b][/color]
  [COLOR=#804040][b]end do[/b][/color]  

  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Rank of the Matrix ='[/color], rank

[COLOR=#6a5acd]9[/color] [COLOR=#804040][b]format[/b][/color]([COLOR=#ff00ff]'s['[/color], i1, [COLOR=#ff00ff]']= '[/color], [COLOR=#008080]f27.20[/color])  
[COLOR=#a020f0]end program[/color] RankOfMatrix

$ gfortran -o happykaren happykaren.f95 -L. liblapack.dll

$ happykaren
 lwmax   =           25
 work(1) =    335.00000000000000     
 lwork =           25
 info =            0
 Singular values were succesfully computed:
s[1]=     35.32704346531138384080
s[2]=     19.99999999999999289457
s[3]=     19.59591794226542305069
s[4]=      0.00000000000000206846
s[5]=      0.00000000000000055535
s[6]=      0.00000000000000000000
s[7]=      0.00000000000000000000
s[8]=      0.00000000000000000000

 epsilon =   4.99999985859034268E-010
 Rank of the Matrix =           3

 

[mikrom]

Now I see: the example program I posted above computes the Rank properly only for non negative singular values.
Use the ABS function.

 

[tesaoghost]

How should I use LAPACK?
Where should I copy liblapack.lib and liblapack.dll files? I get these files but I don't how can I use hem!!
And I am using microsoft visual fortran 2010.

 

[mikrom]

How should I use LAPACK?
...
And I am using microsoft visual fortran 2010.

I don't have knowledge about compiler which you are using. But maybe happykaren100 has found out how to do it.

 

[tesaoghost]

which version of fortran are you using mikrom?
Is your version compatible with windows 7?
I was using microsoft powerstation 4, but this version is not compatible with windows 7, so one of my friends give me microsoft visual fortran 2010 and it is compatible with windows 7. Can I have the version you work with it? Would you please send me its download link?
thanks a lot dear mikrom

 

[mikrom]

which version of fortran are you using mikrom?
Is your version compatible with windows 7?
...
Can I have the version you work with it? Would you please send me its download link?

As you could see above, I compiled the example with gfortran.
It runs on Windows 7.
You can download the compiler for free.

 

[happykaren100]

In fact, I have not solved the problem mentioned yet...
I have tried to use the gfortran, but I still could not solve the problem since the installation was so complicated.
what's more, the gfortran even created more problems in my program as it did not support the library "Portlib"

I am now trying hard to solve the problems in using gfortran...

 

[mikrom]

I compiled the example given above with gfortran, but I thought that you will try to link the Lapack DLL with your compiler.
IMHO, it's possible - you only have to read manual of your compiler: which command line switches you have to use, or where in the IDE you have to make appropriate settings.
But you need first to know, what compiler you are using. Note, that MS Studio isn't a fortran compiler - it's an IDE. Maybe more fortran compilers have interface to run from MS Studio. But I have no experience with this topics because I don't use this IDE.

 

[mikrom]

Hi happykaren100 and tesaoghost,

Here on the LAPACK for Windows page there is an HOWTO: Running LAPACK under Windows.
In the Part 2: Using LAPACK subroutines in a Visual Studio INTEL FORTRAN Project there is an example projekt to download: LAPACK-VS-Example.zip. Have you tried it?