Problem with the rank of the matrix

happykaren100 · Jan 12, 2014

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 · Jan 12, 2014

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 · Jan 12, 2014

Other example - for DGESVD is here:

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

happykaren100 · Jan 15, 2014

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 · Jan 15, 2014

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 · Jan 16, 2014

Hi Happy Karen,

Here is the working example:

Code:

[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

Code:

$ 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 · Jan 16, 2014

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

tesaoghost · Jan 30, 2014

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 · Jan 30, 2014

tesaoghost said:
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 · Feb 2, 2014

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 · Feb 2, 2014

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 · Feb 7, 2014

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 · Feb 8, 2014

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 · Feb 9, 2014

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?

Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

Problem with the rank of the matrix

happykaren100

Technical User

mikrom

Programmer

mikrom

Programmer

happykaren100

Technical User

mikrom

Programmer

mikrom

Programmer

mikrom

Programmer

tesaoghost

Programmer

mikrom

Programmer

tesaoghost

Programmer

mikrom

Programmer

happykaren100

Technical User

mikrom

Programmer

mikrom

Programmer

Similar threads

Part and Inventory Search

Sponsor