Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Highlighted

charlotte_M_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-20-2016
02:53 PM

46 Views

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 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

Accepted Solutions

Highlighted

jimdempseyatthecove

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
11:22 AM

42 Views

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

22 Replies

Highlighted

andrew_4619

Valued Contributor III

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-20-2016
03:53 PM

42 Views

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!

:-)

Highlighted
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.

andrew_4619

Valued Contributor III

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-20-2016
04:02 PM

42 Views

Highlighted

charlotte_M_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
04:11 AM

42 Views

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

Highlighted

andrew_4619

Valued Contributor III

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
05:54 AM

42 Views

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.....

Highlighted

charlotte_M_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
07:32 AM

42 Views

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).

Highlighted

jimdempseyatthecove

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
08:05 AM

42 Views

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

Highlighted

charlotte_M_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
08:39 AM

42 Views

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).

Highlighted

charlotte_M_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
08:53 AM

42 Views

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

Highlighted

charlotte_M_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
10:04 AM

42 Views

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?

Highlighted

jimdempseyatthecove

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
11:22 AM

43 Views

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

Highlighted

jimdempseyatthecove

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
11:27 AM

42 Views

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

Highlighted

jimdempseyatthecove

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
11:41 AM

42 Views

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

Highlighted

charlotte_M_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
11:44 AM

42 Views

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

Highlighted

jimdempseyatthecove

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
12:11 PM

42 Views

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

Highlighted

charlotte_M_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
12:18 PM

42 Views

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.

Highlighted

andrew_4619

Valued Contributor III

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
12:19 PM

42 Views

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..... `

Highlighted

charlotte_M_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
12:30 PM

42 Views

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

Highlighted

charlotte_M_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
12:39 PM

42 Views

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)

Highlighted

charlotte_M_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

09-21-2016
01:43 PM

42 Views

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

For more complete information about compiler optimizations, see our Optimization Notice.