topic Searching online, i think my in Intel® Fortran Compiler

Time splitting issue

charlotte_M_ — Tue, 20 Sep 2016 21:53:06 GMT

Hello Everyone,

I'm trying to solve a problem use time splitting method. For this i have to solve linear systems of equation using tridag (fortran ) routine.

I have two problems

The first one is that the coefficients use in tridag are 2 D array while this routine use 1 D array.

What i did to overcome this is the following (just a example):

do k=1,NT ! Main loop over the time step

!Compute the solution of the first equation

do i=1, N1

do j= 1, N2

A(j) = x(i)*y(j)

B(j)= x(i)*y(j)

C(j)=x(i)*y(j)

End do

Call tridag( A,B,C,F0, F1, N) !F0 is a two dimensional function , F1, a one dimensional function

End do

!Once done i will use F1 as the initial condition for the second equation

Do j =1, N2

Do i=1,N1

E(i)= xx(i)*y(jj)

F(i) =xx(i)*yy(j)

G(i)= xx(i)*yy(i)

End do

Call tridag(E,F,G,F1,F2,n1)

end do

!Next step Write result in the file, for different time

!Finally i update my initial condition

DO i=1,N1

DO j=1,N2

F0(i,j) =F2(i,j)

end do

END DO !ENd of the time loop

****** The second problem is that , i have the same result for different time that i print the solution.

I am doing something wrong?

sorry for the long text

Thank you

do k=1,NT ! Main

andrew_4619 — Tue, 20 Sep 2016 22:53:00 GMT

do  k=1,NT         ! Main loop over the time step
   !Compute  the solution of the first equation 
   do i=1, N1
      do j= 1, N2
         A(j) = x(i)*y(j)
         B(j)= x(i)*y(j)
         C(j)=x(i)*y(j)
      End do 
      Call tridag( A,B,C,F0, F1, N) 
      !F0  is a two dimensional function , F1, a one dimensional function
   End do 

   !Once done i will use F1 as the initial condition for the second equation
   Do j =1, N2
      Do i=1,N1
         E(i)= xx(i)*y(jj)
         F(i) =xx(i)*yy(j)
         G(i)= xx(i)*yy(i)
      End do 
      Call tridag(E,F,G,F1,F2,n1)
   end do 

   !Next step Write result in the file, for different time
   !Finally   i update my initial condition
   DO i=1,N1
      DO j=1,N2
         F0(i,j) =F2(i,j)
      end do
   end do
END DO         !ENd of the time loop

! if you put the code in a Fortran window ( the {code} icon) from the tool bar it is more readable!

:-)

I am slightly confused there

andrew_4619 — Tue, 20 Sep 2016 23:02:34 GMT

I am slightly confused there seem to be several sources (variants) of TRIDAG from a quick google search but all the ones I looked at have 1D array inputs, perhaps you could ID the source of TRIDAG you are using or post the first few lines with the definitions of the items passed in/out.

Hello Andrew_4619,

charlotte_M_ — Wed, 21 Sep 2016 11:11:01 GMT

Hello Andrew_4619,

Thank you for the answer and the fortran code

This is the routine.

SUBROUTINE tridag(a,b,c,r,u,n) 
PARAMETER (nmax=100) 
REAL gam(nmax),a(n),b(n),c(n),u(n),r(n) 
INTEGER j,n 
if (b(1)==0.) pause 'b(1)=0 in tridag' 
bet=b(1) 
u(1)=r(1)/bet 
do j=2,n 
    gam(j)=c(j-1)/bet 
    bet=b(j)-a(j)*gam(j) 
    if (bet==0.) pause 'bet=0 in tridag' 
    u(j)=(r(j)-a(j)*u(j-1))/bet 
end do 
do j=n-1,1,-1 
    u(j)=u(j)-gam(j+1)*u(j+1) 
end do 
END SUBROUTINE tridag

OK I do not know what the

andrew_4619 — Wed, 21 Sep 2016 12:54:40 GMT

OK I do not know what the program is trying to do but you are passing a 2D array to a 1D array so in the first call to TRIDAG you are using only the first N items of the array which is in column order ie F0(1,1), F0(2,1),....F(N,1) or if the size first index of F0 is <N you will wrap into column 2.

I am guessing you want to pass a slice of F0 e.g. F0(:,i) (passes a 1D array that is the column number i) or maybe F0(j,:) for row number J. Note that the second case is not effecient as the slice represents non-adjacent memory locations so the compiler will create a temporary array to hold this data.....

For the first equation What

charlotte_M_ — Wed, 21 Sep 2016 14:32:28 GMT

For the first equation What i want is for a given value of i, solve the equation for all the value of j.

The process has to be repeat for all the value of i ( loop over i).

The initial distribution is F0(i,j).

Your comment in your code

jimdempseyatthecove — Wed, 21 Sep 2016 15:05:52 GMT

Your comment in your code posting state that F0 and F1 are functions, not arrays. The subroutine tridag specifies that dummy arguments r and u are REAL arrays.

If your comment is wrong, or should I say partly wrong, in that F0 and F1 are arrays, with F0 being a 2D array, then Andrew's comment about what are you intending to pass into tridag as r? IOW, is it a row in F0 or a column in F0?

Your code having the outer DO k=1,NT would seem to indicate that F0 may contain different time steps as either by row or by column. We do not have information as to how F0 is arranged.

Jim Dempsey

charlotte_M_ — Wed, 21 Sep 2016 15:39:43 GMT

Sorry my comment is wrong. F0(i,j) and F1 are both real arrays. F0 is 2 D array.

F0 is r in tridag.

Let's say F0(i,j) = exp(- i**2)*exp(-j).

Searching online, i think my

charlotte_M_ — Wed, 21 Sep 2016 15:53:33 GMT

Searching online, i think my problem is how do i pass this 2 D array to tridag routine?

And once i have the result from tridag, how can i have back the array.

Thank you so much

In all the case i try to pass

charlotte_M_ — Wed, 21 Sep 2016 17:04:09 GMT

In all the case i try to pass in tridag the index of the second loop .. It means a column in F0.

However for the second equation, the column becomes the row and the row the column.

Little mistake, the first call is CALL tridag( A,B,C,F0,F1,N2)

F0 is allocate array of dimension (N1, N2)..

So yes i want to pass a slice of F0

For the first equation i want to pass F0(i,:) and the second F0(:,j)

Any suggestion?

Pass F0(i,:) and F0(:,j).

jimdempseyatthecove — Wed, 21 Sep 2016 18:22:18 GMT

Pass F0(i,:) and F0(:,j).

Note, in the case of F0(i,:), the compiler will generate a temporary (because the slice has a stride of N1, size of N2).
For the F0(:,J), the compiler should use the data in place (because the slice has a stride of 1, size of N1).

Jim Dempsey

By the way. Consider using:

jimdempseyatthecove — Wed, 21 Sep 2016 18:27:00 GMT

By the way. Consider using:

SUBROUTINE tridag(a,b,c,r,u,n) 
REAL gam(n),a(n),b(n),c(n),u(n),r(n) ! have gam assume the size of the other arrays (n)

This will prevent a buffer overflow should n (N1 or N2 from caller) exceed 100.

If n is very large (say over 1,000,000) consider making gam allocatable:

SUBROUTINE tridag(a,b,c,r,u,n) 
REAL a(n),b(n),c(n),u(n),r(n)
REAL, allocatable :: gam(:)
allocate(gam(n))
...
deallocate(gam) ! optional on newer specifications

Jim Dempsey

Once you get the code working

jimdempseyatthecove — Wed, 21 Sep 2016 18:41:10 GMT

Once you get the code working, and if the entire process takes longer to execute than you would like, readers here can show you how to get better performance.

If N1 and N2 get reasonably large, then parallelization and/or vectorizaton improvements may be available.
And, even if parallelization is not warranted, two threads can have one thread computing time steps and the other thread performing the file writes. IOW runtime goes from:

(computation time) + (File Write time)

to:

max((computation time), (File Write time)) + (a tad for threading overhead)

Jim Dempsey

Thank you so much Jim

charlotte_M_ — Wed, 21 Sep 2016 18:44:39 GMT

Thank you so much Jim Dempsey for your answer.

For the first equation i passed F0(i,:). Everything was ok.

The solution given by TRIDAG is F1.

Now i want to use F1 as F0 for the next call (second equation). F1 becomes my initial array (IOW r)

But F1 is not a 2D array, how can i pass F1(j,:)

I made the modification for the gam ..Thank you so much

In looking at your original

jimdempseyatthecove — Wed, 21 Sep 2016 19:11:14 GMT

In looking at your original sketch code (#1), you are not using the F1 results after the call to TRIDAG and its end do. IOW you are overwriting the results. Therefor, I assume you intended to have made F1 2D as well, then process the row/col of F1 array in the second nested loop.

Jim Dempsey

Quote:jimdempseyatthecove

charlotte_M_ — Wed, 21 Sep 2016 19:18:28 GMT

jimdempseyatthecove wrote:

In looking at your original sketch code (#1), you are not using the F1 results after the call to TRIDAG and its end do. IOW you are overwriting the results. Therefor, I assume you intended to have made F1 2D as well, then process the row/col of F1 array in the second nested loop.

Jim Dempsey

Yes exactly, after making F1 2D array , i pass F1(:,j) for the second call.

well in the first case you

andrew_4619 — Wed, 21 Sep 2016 19:19:00 GMT

well in the first case you have Call tridag( A,B,C,F0(i,:), F1, N2)

but your second points confuses me, is it that at this first call you need the values of F1 for each i to use later?

If so make F1 2D to the same dims as F0 and then have as Call tridag( A,B,C,F0(i,:), F1(i,:), N2)

then later you would have: Call tridag(E,F,G,F1(:,J),F2,n1)

I am guessing as I have made no attempt to understand teh nature of your calculation.....

F1 is the solution of

charlotte_M_ — Wed, 21 Sep 2016 19:30:58 GMT

F1 is the solution of TRIDAG, no need to know it.

When i call TRIDAG for the first time , the unknown is F1, so TRIDAG will give me F1

This is what i did :

charlotte_M_ — Wed, 21 Sep 2016 19:39:47 GMT

This is what i did :

First call : Call TRIDAG ( A,B,C,F0(i,:), F1(i,:) , N2)

Then for the second

call TRIDAG (E,F,G, F1(:,j), F2(:,j) , N1)

charlotte_M_ — Wed, 21 Sep 2016 20:43:37 GMT

Thank you so much Jim Dempsey and andrew_4619.

Everything is working now. And i have a better understanding.

The next step will be parallelization..

So see you soon