NaN

[ferry2]

I just finished my application but for this subroutine

subroutine back_substitution(s, c)
implicit none

double precision, dimension,, intent(out) :: s
double precision, dimension), intent(out) :: c

real :: w
integer :: i, j, n

n = size(c)


do i = n, 1, -1
w = c(i)
do j = i + 1, n
w = w - s(i, j)*c(j)
end do
c(i) = w/s(i, i)
end do

end subroutine back_substitution

I'm receiving NaN for c(i).

Why is that happen?

 

[mikrom]

C has size of N.
But, when you have in 1.loop i=N, then in 2. loop for j=N+1 you try to compute C(N+1) which is undefined.

 

[FJacq]

I think that Mikrom is wrong because the second loop goes from i+1 to n. So the loop is never executed when i==n.

But I remark a mistake anyway : you have declared the arrays "c" and "s" with INTENT(out). As they are not initialized in your subroutine before being used (see the instruction w=c(i) for instance), then they could contains a lot of undefined components (the compiler may fill them with anything).

I assume that these arrays are normally initialized outside the subroutine. So the array "c" must be declared INTENT(inout) and the array "s", because not modified, INTENT(in).

François Jacq

 

[ferry2]

I already tryed this @FJacq. I still can't understand why I'm geting this result.

 

[salgerman]

a couple of print statements could go along way...try to see what the content of s() and c() as you enter the subroutine and before you do any calculations..

 

[FJacq]

You should print out the matrix "s" and the vector "c" at the entrance of the subroutine. You have perhaps s(i,i) almost equal to zero for a given i. This is often the case in a Gauss method when the matrix is almost singular (your programming looks like the final step of the Gauss algorithm).

Don't forget to correct the intent anyway...

François Jacq

 

[ferry2]

Hi, @FJacq. You correctly suggests that this is part of gaussian ellimination. I do what you told me. I correct the intent of c to (in out). I print the matrix "s" and the vector "c" at the entrance of the sub routine and the result is only zeros.

 

[salgerman]

Then, you have to take one step back and see why 's' and 'c' are zeros! I presume that if you print 's' and 'c' (or the variables that you are using in the call to back_substitution) you are also getting zeroes?

...keep backtracking your steps until you why you have zeros...

 

[FJacq]

Did you change also the INTENT for the matrix S (I suggested INTENT(IN)) ?

You indicate that you have changed only the INTENT for C : this is not enough of course ! With INTENT(out), your matrix could be filled with 0 automatically ...

François Jacq

 

[ferry2]

Yes, I change the intent for s to intent(in). Sorry. I forgot to mention that.

 

[FJacq]

So print s and c just BEFORE the call statement to the routine back_substitution...

This is a detective work... You are Sherlock Holmes. Try to locate the origin of the trouble.

François Jacq

 

[FJacq]

Another question : is your routine back_substitution inserted in a module ?

This is important because you have declared the argument arrays with the dimensions , and ). These declarations mean that the array sizes are passed automagically and may be get with the intrinsic function SIZE.

But this mode of argument passing must be known by the compiler when analyzing the procedure calling your subroutine. The easiest way to indicate this mode consists in having a USE statement followed by the name of the module containing your subroutine. Without such use statement, the default mode of argument passing is the Fortran-77 one which is not the mode ,.

Notice that if the calling procedure and the called subroutine are put inside the same module, then the compiler knows already everything and the USE statement is useless and even forbidden : a module cannot use itself !.

François Jacq

 

[ferry2]

Hi, @FJacq.
Yes the "back_substitution" subroutine is in module file which I call in main program with USE statement.
i think the problem with NaN value is caused by the fact that the matrix have too much zeros and some where may be occurs dividing by zero.
I have to implement some other alghorithm form solving linear system of equations.

 

[FJacq]

If your matrix is singular, then there is no algorithm able to solve "exactly" your linear system because it has no solution at all (or an infinite number of solutions)...

Only a minimization algorithm could give the "nearest solution".

Notice that the Gauss algorithm *IS* the main direct algorithm to solve a general linear system. For specific matrices, other algorithms could apply. For instance, with symmetrical definite positive matrices, The Cholesky algorithm is preferred. For very large sparse matrices, iterative algorithms may be more suitable (for instance conjugate gradient with Cholesky preconditioning if the matrix is also symmetrical).


François Jacq