using functions for tensor math

[titous09]

hello,

i'm trying to creat a function that will take two (3x3) matrices as an input and output a (3x3x3x3) tensor. i have already created a code, but when i try and compile i get the following error: Error: The number of subscripts is incorrect.

Here is my source code:
program tensors
IMPLICIT NONE


REAL, DIMENSION(3,3) :: s
REAL, DIMENSION(3,3,3,3) :: J, doubledot


call kronecker(s)
call print2(s)

J=doubledot(s,s)

END



SUBROUTINE kronecker(a)

REAL, DIMENSION(3,3) :: a
INTEGER :: w, x

DO w=1,3
DO x=1,3
IF (w==x) THEN
a(w,x)=1
ELSE
a(w,x)=0
END IF
END DO
END DO

END



FUNCTION doubledot(a,b)

REAL, DIMENSION(3,3) :: a,b
REAL, DIMENSION(3,3,3,3) :: doubledot
INTEGER :: w, x, y, z


DO w=1,3
DO x=1,3
DO y=1,3
DO z=1,3
doubledot(w,x,y,z)=a(w,x)*b(y,z)
END DO
END DO
END DO
END DO


END


SUBROUTINE print2(a)

REAL, DIMENSION(3,3) :: a
INTEGER :: w

101 FORMAT (3F6.2,1X)
DO w=1,3
PRINT 101, a(w,1),a(w,2),a(w,3)
END DO

END


SUBROUTINE print4(a)

REAL, DIMENSION(3,3,3,3) :: a
INTEGER :: w, x, y, z


102 FORMAT (F4.2,1X,A,1X,4F2.0)
DO w=1,3
DO x=1,3
DO y=1,3
DO z=1,3
PRINT 102, a(w,x,y,z),"at",w,x,y,z
END DO
END DO
END DO
END DO

END


i feel like this is just a simple error that i can't get around. thanks for the help!

 

[titous09]

well i found a way around this. you have to use the "interface" structure. do something like this:
interface doubledot
function doubledot(a,b)
real, dimension(3,3) :: a, b
real, dimension(3,3,3,3) :: doubledot
end function doubledot
end interface

thanks for the help though!

 

[gummibaer]

Well,

I am not smart enough to understand this thing with the interface, but the bug you have is, that in your main program you define doubledot to be a variable of rank 4 (=number of subscripts)

REAL, DIMENSION(3,3,3,3) :: J, doubledot

and a few lines further down you use doubledot with two subscripts only.

J=doubledot(s,s)

In this you use the entity J of rank 4 to be equal to doubledot with the wrong number of subscripts - which will not work.

Next: In main you define doubledot to be a variable of rank 4. But later you define a function under the same name. I guess, your function doubledot will never be executed because doubledot is declared a variable. Your linker should have issued a warning of this multiple definition of different entities under the same name.

I gues, for what you want to do, write

real,dimension(3,3,3,3)::J

...

call doubledot (J,s,s)

...

END


Soubroutine doubledot (J,a,b)
REAL, DIMENSION(3,3) :: a,b
REAL, DIMENSION(3,3,3,3) :: J
INTEGER :: w, x, y, z


DO w=1,3
   DO x=1,3
       DO y=1,3
           DO z=1,3
           J(w,x,y,z)=a(w,x)*b(y,z)
           END DO
       END DO
   END DO
END DO


END