couple of observations, Fortran needs matrix diagonal() function, and determinant()

For fun, I spend over 3 hrs adding Fortran to one entry
of my HOWTO cheat sheet page. This entry is on how to compute
the adjugate matrix of a square matrix.

The adjugate matrix is defined here

http://en.wikipedia.org/wiki/Adjugate_matrix

It is basically the transpose of the cofactor matrix. If you
divide it by the determinant of the matrix, you get the
inverse of the matrix.

I thought this will take me may be 10-20 minutes in Fortran,
but it ended up taking more than 3 hrs. The first reason is
that Fortran does not have determinant function build-in,
and also, it does not have a function to read the diagonal
of a matrix. So, ended up linking to lapack DGETRF() to
do LU decomposition in order to obtain the determinant
from the PLU result of the decomposition.

But all went well in the end. I get the same results as
shown on the wiki, and the same as Matlab and Mathematica,
tested on number of matrices and learned a bit of Fortran.

But my question is: Why is there no intrinsic for finding
the determinant for rank-2 matrix? I find this strange
given Fortran is for numerical computation and given
that Lapack does not have a build-in function either
(one needs to use LU and more processing afterwords)

And also there is no function to obtain the diagonal of
rank 2 array? I see Fortran now has intrinsic for matrix
multiply, transpose. I think a diag() function would be
nice to have, and also a det() function.

I found Fortran array/matrix operation to be flexible
and easy to work with. I just think it needs few more
intrinsic functions such as the above missing ones to make
it 'easier' to do common things. For example, I did this
same operation in Mathematica in only 3 lines, a solution in
Matlab is the same.

Here is the Fortran function I wrote. If you find improvements
please feel free to let me know and I'll correct

-------------------------------------
!-- Find the Adjugate of a square matrix
!-- Version 6/22/2012
!-- gfortran 4.6.3
program t44
implicit none

integer, parameter :: n=3
integer :: INFO,i,j,ii,IPIV(n-1),number_row_exchange
real (kind=kind(0.0d0)) :: A(n,n),B(n-1,n-1),C(n,n)
logical :: M(n,n)
A(1,:) = [1,2,3];
A(2,:) = [4,5,6];
A(3,:) = [7,8,9];

!A(1,:) = [-3,2,-5];
!A(2,:) = [-1,0,-2];
!A(3,:) = [3,-4,1];
DO i=1,n
DO j=1,n
!build the mask to extract submatrices
M = .true.
M(:,j) = .false.
M(i,:) = .false.
B = reshape(pack(A, M),[n-1,n-1])
!-- LU decomposition in order to obtain determinant
CALL DGETRF( n-1, n-1, B, n-1, IPIV, INFO )
!-- determinant of U is the product of diagonal
C(i,j)=1.0d0
DO ii=1,n-1
C(i,j)=C(i,j)*B(ii,ii)
END DO
!-- Adjust sign based on number of row exchanges
number_row_exchange = 0
DO ii=1,n-1
IF (IPIV(ii) .NE. ii) THEN
number_row_exchange = number_row_exchange +1
END IF
END DO
IF (MOD(number_row_exchange,2).EQ.1) THEN
C(i,j)=-C(i,j)
END IF

!-- Apply the (-1)^(i+j) logic
IF (MOD(i+j,2).EQ.1) THEN
C(i,j)=-C(i,j)
END IF
END DO
END DO

!-- Transpose the matrix of cofactor to obtain the adjugate matrix
C = TRANSPOSE(C)

PRINT *,'C='
write(*,'(3F15.4)') ( (C(i,j),j=1,n), i=1,n)

end program t44
-----------------------------

compile/link/run (linux)

export LD_LIBRARY_PATH=/usr/lib/atlas-base/:/usr/lib/atlas-base/atlas/
gfortran -fcheck=all -Wall t44.f90 -L/usr/lib/atlas-base/atlas/ -lblas -llapack

./a.out
C=
-3.0000 6.0000 -3.0000
6.0000 -12.0000 6.0000
-3.0000 6.0000 -3.0000

fyi, using Mathematica, has build-in Cofactor function:
-------------------------
<<"Combinatorica`"
mat = {{1,2,3},{4,5,6},{7,8,9}};
cof = Table[Cofactor[mat,{i,j}],{i,3},{j,3}];
adjugate = Transpose[cof]
------------------------

Out[22]= { {-3, 6, -3},
{ 6, -12, 6},
{-3, 6, -3}}

--Nasser
.

Relevant Pages

  • Re: Strange Fortran version
    ... "ENTRY lable" was common in some Fortran IV subroutines and probably ... Fortran compilers (many were written as a higher exercise in M.SC tasks ...
    (comp.lang.fortran)
  • Re: ENTRY statement
    ... < subroutine at its two different entry points. ... Where they are in memory, ... became the Fortran 66 standard. ... one was allowed to access dummy arguments from entry points ...
    (comp.lang.fortran)
  • Re: ENTRY statement
    ... produce libraries that contain a Fortran77 style alias (for example I ...    end interface ... These all have some advantages and disadvantages compared to entry ... For procedures to be used by Fortran 90 and higher compilers these ...
    (comp.lang.fortran)
  • Re: Problem with Return Value from a Function
    ... it presumably came from Fortran. ... EQUIVALENCEs the return values of all function ENTRY points, ... Compilers I know of generate an epilog routine for each ... there are no ENTRY statements it should be easy to do at compile time. ...
    (comp.lang.pl1)
  • Re: ENTRY statement
    ... I thought an ENTRY statement ... As others have noted, x, being a dummy argument of another entry ... Fortran doesn't specify call by reference, but does allow for it, ... the subprogram need not appear in the argument list of the ...
    (comp.lang.fortran)