Hi,

I'm a beginner in NR and also a relatively beginner in C++.
I'm basically trying to use LU decompostion to solve a set of linear equations. (Using ludcmp and lubksb)
But I'm just not sure how to call functions and send appropriate values to the right variables.
Like I said, I'm relatively a newbie to both NR and C++.

This is how my code is laid out:

int main()
{
/*User inputs certain values then it carries out a few calculations to obtain the matrix that is sized dim x dim. Matrix is called 'a'. Also creates the right hand side matrix (a single column) called 'b'.*/
return 0;
}

void ludcmp(float **a, int dim, int *indx, float *d)
{
/*I guess **a refers to the matrix I created in the main function, and dim is the size of the matrix, *indx is a vector that is needed to carry out the inversion process, same thing with d. At this point, I'm not sure if my matrix 'a' gets sent to the ludcmp. Same concern for 'dim'*/
}

void lubksb(float **a, int dim, int *indx, float b[])
{
/*In here, **a is NOT supposed to be the matrix I create in the main function, it is supposed to be a matrix modified in ludcmp. I'm not sure if things work properly here. 'dim' is supposed to be the size of the matrix as I mentioned above. b[] is the right hand side matrix created in the main function. I'm not sure if that gets sent in here as well. Finally, the output x (solution) gets written over the matrix b[].. if I'm not mistaken. So I need to output this as well some how.*/
}

int main2()
{
/*I call out the functions here by doing something like,

float **a,*b,d;
int dim,*indx;
...
ludcmp(a,dim,indx,&d);
lubksb(a,dim,indx,b);

but again, I'm not sure if variables go into right places and gets called by appropriate items.*/
}



So there you have it, I hope the explanation was good enough.
If someone helps me with this, I'll be forever thankful to you.

PK

First you say

... in C++...
Then you show us

void ludcmp(float **a, int dim, int *indx, float *d)
and

void lubksb(float **a, int dim, int *indx, float b[])

Those are the function prototypes for the C version. Are you doing C or C++? There is a difference in how a main() program declares and populates its variables and calls the functions. (Big difference.)

Regards,

Dave

I'm using C++.
I didn't think there was a huge difference between them.
I refer to http://www.cplusplus.com/doc/tutorial/ mostly, and they had the same structure for function prototypes?

Could you show me the proper way for C++?

Let me give you what I've got so far:

//HEADER FILES
#include <iostream>
#include <math.h>
#include <cmath>
#include <fstream>
#include "nrutil.h"
#include "nrutil.c"
#define TINY 1.0e-20;
using namespace std;

//LU Decomposition - straight off the book with a minor change of 'n' becoming 'dim' (matrix size)
void ludcmp(float **a, int dim, int *indx, float *d)
{
int i,imax,j,k;
float big,dum,sum,temp;
float *vv;

vv=vector(1,dim);
*d=1.0;

for (i=1;i<=dim;i++)
{
big=0.0;
for (j=1;j<=dim;j++)
if ((temp=fabs(a[i][j])) > big) big=temp;
if (big == 0.0) nrerror("Singular matrix in routine ludcmp");
vv[i]=1.0/big;
}

for (j=1;j<=dim;j++)
{
for (i=1;i<j;i++)
{
sum=a[i][j];
for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.0;
for (i=j;i<=dim;i++)
{
sum=a[i][j];
for (k=1;k<j;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (dum=vv[i]*fabs(sum)) >= big)
{
big=dum;
imax=i;
}
}
if (j != imax)
{
for (k=1;k<=dim;k++)
{
dum=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=dum;
}
*d = -(*d);
vv[imax]=vv[j];
}
indx[j]=imax;
if (a[j][j] == 0.0) a[j][j]=TINY;
if (j != dim)
{
dum=1.0/(a[j][j]);
for (i=j+1;i<=dim;i++) a[i][j] *= dum;
}
}
free_vector(vv,1,dim);
}

void lubksb(float **a, int dim, int *indx, float b[])
{
int i,ii=0,ip,j;
float sum;
for (i=1;i<=dim;i++)
{
ip=indx[i];
sum=b[ip];
b[ip]=b[i];
if (ii)
for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
else if (sum) ii=i;
b[i]=sum;
}
for (i=dim;i>=1;i--)
{
sum=b[i];
for (j=i+1;j<=dim;j++) sum -= a[i][j]*b[j];
b[i]=sum/a[i][i];

for (i=1; i<=dim; i++)
{
cout<<b[i]<<"\n";
}
}
}

//Main function
int main ()
{
float **a,d,*b;
int *indx;
int i,j,k,n,m,r,f,s,g,dim;

...
...
...

ludcmp(a,dim,indx,&d);
lubksb(a,dim,indx,b);

system ("pause");
return 0;
}


Thanks.

PK

I'm using C++.
I didn't think there was a huge difference between them.
There are differences that make some C code not compatible with C++. (Some compilers are more picky than others.) It may be possible to compile and use the particular C functions that you use in a C++ program, and I show an example. If your C++ compiler complains and you can't resolve the problem, maybe you could try changing your main() function file to C (using printf, fscanf, etc., instead of cin, cout, etc.)

As for your functions, there is a bug in the code that you added to lubksb to print out the solution vector:


void lubksb(float **a, int dim, int *indx, float b[])
{
.
.
.
for (i = dim; i >= 1; i--) {
sum = b[i];
for (j = i + 1; j <= dim; j++)
sum -= a[i][j] * b[j];
b[i] = sum / a[i][i];

for (i = 1; i <= dim; i++) { /*<--- This is the wrong place */
cout << b[i] << "\n";
}
}
}

I made the indentation consistent with the way that I like to see it. (Of course C and C++ compilers ignore indentation, but it helps me keep things straight.)

The print loop has to be outside the loop that is calculating the elements of b. I'll put it in the right place (with a little more narrative). In fact, I'll just comment it out, since I don't think that this function is the proper place for output except for debugging. (If you want to see it here, then change "#if 0" to "#if 1")


for (i = dim; i >= 1; i--) {
sum = b[i];
for (j = i + 1; j <= dim; j++)
sum -= a[i][j] * b[j];
b[i] = sum / a[i][i];
}
#if 0
/* for debugging change the above from "#if 0" to "#if 1" */
for (i = 1; i <= dim; i++) {
cout << "b[" << i << "] = " << b[i] << endl;
}
cout << endl;
cout << "Leaving lubksp" << endl;
#endif



}

Now, back to the main part of your question. How to get stuff into the arrays and use the functions.

Here's my test program in C++, using nrutil.c and your functions:

/* The input file contains the following:
* First line is two integers. The first is
* the size of the system. The second is
* the number of right-hand side vectors.
*
* The variables for these are n and m.
*
* Everything else on the first line is ignored.
*
* Then there are n*n floats, separated by
* whitespace. There can be any number of
* items on each line as long as there are
* a total of n*n values.
*
* Then there are m right-hand side vectors
* Each one has n floats, separated by whitespace.
*
*
* A sample input file, indata.txt, might look like:

3 2 #Size of system and number of right hand sides
1 2 3
2 3 4
1 3 4
6 9 8
14 20 19

Output will be something like

File indata.txt: 3x3 matrix

1.000000e+00 2.000000e+00 3.000000e+00
2.000000e+00 3.000000e+00 4.000000e+00
1.000000e+00 3.000000e+00 4.000000e+00

=============================================
Solution:
+1.000000e+00 +1.000000e+00 +1.000000e+00
Right-hand side:
+6.000000e+00 +9.000000e+00 +8.000000e+00
Matrix * solution:
+6.000000e+00 +9.000000e+00 +8.000000e+00
=============================================
Solution:
+1.000000e+00 +2.000000e+00 +3.000000e+00
Right-hand side:
+1.400000e+01 +2.000000e+01 +1.900000e+01
Matrix * solution:
+1.400000e+01 +2.000000e+01 +1.900000e+01
=============================================

*/

#include <iostream>
#include <iomanip>
#include <fstream>
#include <cstdlib>
#include <cmath>
#define TINY 1.0e-20;
#define NRANSI

#include "nrutil.h"
#include "nrutil.c"
using namespace std;

/*ludcmp and lubskp functions can go here */

int main(int argc, char **argv)
{
int i, j, k, m, n, *indx;
int num;
float p, *x, **a, **aa, *b;
float sum;
ifstream infile;

if (argc < 2) {
cout << "Usage: " << argv[0] << " filename" << endl;
exit(EXIT_FAILURE);
}
infile.open(argv[1]);
if (!infile) {
cout << "Cannot open file " << argv[1] << " for reading." << endl;
exit(EXIT_FAILURE);
}
cout << "File " << argv[1];
if (!(infile>>n>>m)){
cout << "\nFirst line must have size of system and number of rhs\n";
exit(EXIT_FAILURE);
}
if (m < 1) {
cout << "\nMust have at least one right-hand side vector\n";
exit(EXIT_FAILURE);
}
/* make an "empty" loop to ignore the rest of the first line */
while (infile.get() != '\n')
;

indx = ivector(1, n); /* holds pivot info from ludcmp for lubskp*/
a = matrix(1, n, 1, n); /* this holds the original matrix */
aa = matrix(1, n, 1, n); /* a copy of the matrix to send to ludcmp */
b = vector(1, n); /* the original right-hand side */
x = vector(1, n); /* Copy rhs to x and send to lubksp */

cout << ": " << n << "x" << n << " matrix" << endl << endl;

cout << scientific << showpos;

for (i = 1; i <= n; i++) {
for (j = 1; j <= n; j++) {
if (!(infile >> a[i][j])) {
cout << "Couldn't read a[" << i << "][" << j << "]" << endl;
exit(EXIT_FAILURE);
}
aa[i][j] = a[i][j]; /* save a, work with aa */
cout << setw(15) << a[i][j];
}
cout << endl;
}

ludcmp(aa, n, indx, &p); /* perform the decomposition */

for (num = 1; num <= m; num++) { /* loop for the rhs vectors */
cout << "\n=============================================\n";
/* Solve equations for each right-hand vector */
for (k = 1; k <= n; k++) {
if (!(infile >> b[k])) {
cout <<"Couldn't read b[" << k << "]" << endl;;
exit(EXIT_FAILURE);
}
x[k] = b[k]; /* save b, work with x */
}

lubksb(aa, n, indx, x); /* back substitution for this rhs */

cout << "Solution:" << endl;
for (i = 1; i <= n; i++) {
cout << setw(15) << x[i];
}
cout << endl << endl;

cout << "Right-hand side:" << endl;
for (i = 1; i <= n; i++) {
cout << setw(15) << b[i];
}
cout << endl;

cout << "Matrix * solution:" << endl;
for (i = 1; i <= n; i++) {
sum = 0.0;
for (j = 1; j <= n; j++) {
sum += a[i][j] * x[j];
}
cout << setw(15) << sum;
}
}
cout << "\n=============================================\n\ n";

infile.close();

free_vector(b, 1, n);
free_matrix(a, 1, n, 1, n);
free_matrix(aa, 1, n, 1, n);
free_vector(x, 1, n);
free_ivector(indx, 1, n);

return 0;
}



This gets the file name from the user command line. So, if the name of this file is "test.c" and you compile it to be "test.exe" and your input data is a file named "indata.txt" you enter the following command line:

test indata.txt

Here's little test file with a 3x3 matrix and two right-hand side vectors:


3 2 #Size of system and number of right hand sides
1 2 3
2 3 4
1 3 4
6 9 8
14 20 19


Output:

E:\nr2\test>test indata.txt
File indata.txt: 3x3 matrix

1.000000e+00 2.000000e+00 3.000000e+00
2.000000e+00 3.000000e+00 4.000000e+00
1.000000e+00 3.000000e+00 4.000000e+00

=============================================
Solution:
+1.000000e+00 +1.000000e+00 +1.000000e+00

Right-hand side:
+6.000000e+00 +9.000000e+00 +8.000000e+00
Matrix * solution:
+6.000000e+00 +9.000000e+00 +8.000000e+00
=============================================
Solution:
+1.000000e+00 +2.000000e+00 +3.000000e+00

Right-hand side:
+1.400000e+01 +2.000000e+01 +1.900000e+01
Matrix * solution:
+1.400000e+01 +2.000000e+01 +1.900000e+01


Regards,

Dave

I appriciate your help big time Dave.
When my master's gets published next year, I won't forget to mention your name in my 'thanks to' list. ;)

Again, thanks a lot.

PK

Dave,

Thanks to you, I've almost got this baby in bed, however I'm a bit confused at the end.

But for some reason, it doesn't like the first 'if' statement which is:

if (argc < 2) {
cout << "Usage: " << argv[0] << " filename" << endl;
exit(EXIT_FAILURE);
}

and it exits.
What is argc?
I even replaced " filename" to my txt file name, but no luck.

But for some reason, it doesn't like the first 'if' statement
What do you mean it doesn't like it? The purpose of the if statement is to tell you that you must supply a file name on the command line.

So, as I said before:


This gets the file name from the user command line. So, if the name of this file is "test.c" and you compile it to be "test.exe" and your input data is a file named "indata.txt" you enter the following command line:

test indata.txt


If you don't supply a file name, then the main function parameter argc has a value of 1, and the message gets printed

Usage: test filename

In other words, it tells you that you didn't give it a file name to work with.

The idea is that you can run the program with different input files by invoking the program with the name of whatever file you want to use.

The first parameter of the main() function lets the program know how many items were on the command line. The first argument is the name of the program and other arguments consist of whatever else the user entered on the command line.

If the program is going to be run with the same file every time, instead of using argv[1] as the file name, just hard-code it into the program:

int main()
{
int i, j, k, m, n, *indx;
int num;
float p, *x, **a, **aa, *b;
float sum;
ifstream infile;
char * inname = "indata.txt";

infile.open(inname);
if (!infile) {
cout << "Cannot open file " << inname << " for reading." << endl;
exit(EXIT_FAILURE);
}
.
.
.


Now, you have to recompile the program if you want to use another file name. (Which is why I showed it the other way.)

Regards,

Dave

Joly It is finally working!

I can't thank you enough Dave!

Thanks again!

Dave,

It seems like an issue comes up everytime I think the case is closed.
Let me explain what the problem is.

There's actually two parts to my project (two separte exe files); part1.exe creates the matrix that needs to be solved, and part2.exe solves the matrix - this is the part you helped me out with.

Now, at most of the time I'm dealing with BIG matrices, something like 500 x 500 or greater.
However, part1.exe - which is the code that produces the matrix and writes it into txt file - crashes when I try to go beyond let say 300 by 300.
At work, which I've got 3 gb ram, can handle a bit more, up to something like 500 x 500.
At home, I've got a 1.5 gig and it struggles to get passed 300 by 300.

Now I know that this issue is related to the performance of my computer - but is there a way to improve this size by playing around with the code?

Also, what is the max size of n that would be acceptable for the part2.exe - which is the part you helped out with? (in other words, have you tested yours? and if so, what is the max n?)
At work, it works up to 500 as I mentioned, but I wasn't able to test any further because simply my part1 code won't produce the matrix beyond 500.

Thanks always.

part1.exe creates the matrix that needs to be solved,
Now, at most of the time I'm dealing with BIG matrices, something like 500 x 500 or greater.
However, part1.exe - which is the code that produces the matrix and writes it into txt file - crashes when I try to go beyond let say 300 by 300.

Let's see...300x300 floats (assuming 4 bytes per float) is 360,000 bytes. 500x500 floats comes out to 1,000,000 bytes. That is not an outrageous amount of memory. (I won't try to guess at how long it's going to take to solve 500 equations and 500 unknowns---I mean, I could try to guess, but I won't)


At work, which I've got 3 gb ram, can handle a bit more, up to something like 500 x 500.
At home, I've got a 1.5 gig and it struggles to get passed 300 by 300.
Now, the amount of memory available to your C++ program is not directly related to the amount of physical memory in your system. Get out of that mindset. It simply isn't.

Different compilers for different operating systems have different limits, but here's the drill:

The nr stuff (using functions like vector() and matrix() from nrutil.c) uses dynamic memory allocation. The maximum amount of memory that my GNU compilers can get from my 32-bit Linux systems is something like three Gigabytes. My Windows XP compilers (Microsoft, Borland, GNU) can get about 2 Gigabytes worth. I hate to repeat myself, but this has nothing to do with the amount of physical memory on the systems. Also note that this is several orders of magnitude more than you are having problems with.

Now, here's the scoop with C compilers on systems that we see these days: If you declare arrays inside main() or other functions (you know, statements like "float x[500][500];") this gets memory from the program stack. Default settings on my Windows-based compilers crap out at about a Megabyte. (I think that's what I remember; I haven't actually tested lately.) My GNU/LInux systems give a bit more, but the point is that if I am going to require more than a few hundred K bytes of memory in an application, I don't depend on the stack.

One problem is that the program doesn't always tell you what's wrong. It can't tell you at compile time, because it's really a linker limitation, not a compiler roadblock. And at run time, it just quits (maybe with a seg fault or maybe it just silently crashes, or...) They don't always tell you these things in programming classes. Welcome to reality.

Bottom line: don't use the stack for arrays that need hundreds of thousands of bytes.

Some ways around the problem:
1. Use dynamic allocation (C++ new operator or C malloc() function). Maybe even use nr vector() and matrix() function. They do the dynamic allocation for you. And I know for sure that you can get matrices with millions of bytes from them. For sure.

or

2. Use C++ vectors instead of arrays. (And vectors of vectors instead of 2-D arrays) More elegant; more robust for C++ objects; not necessarily more efficient for simple array operations.

I saved the easiest for last:

3. For simple programs (not multithreaded or other complicated multi-tasking stuff with re-entrant signal handlers, etc.) just declare the arrays to be static.

That is, instead of:


int main()
{
.
.
.
float x[500][500];
.
.
.


Just try

int main()
{
.
.
.
static float x[500][500];
.
.
.
}


Now if this doesn't help (or if it doesn't make sense), then how about showing us a little of your code that gives the problem(s)? Tell us what compiler and operating system you are using. (Sometimes it makes a difference to people who want to help.)

Regards,

Dave

I've tried the 'static float' bit, but an error apeared as it states that the float is not constant - my matrix size is not fixed.

Thus I seek your help, and I'll present my code - part1 that is.
If I briefly explain what is going on:

1.User first inputs the number of rows and columns within the flow field (n-rows by m-columns, thus 3 x 4 field gives 12 nodes).

2.Inlet pressure is entered then outlet pressure is entered - take 3 x 4 field for example, the inlet pressures are at the node 1, 5 and 9 - outlet pressures are at the node 4, 8 and 12.
The 'middle' node pressures are unknowns.

3. A set up pressure-related relation developes the equations to solve for the 'middle' unknowns nodes. Thus we have 6 unknowns and 6 equations for this particular case.

4. The equation is written in matrix form to a textfile.

This is the purpose of this code, and the part2 is for solving.
The trouble lies within part1, thus I'll attach part1.
My indentation is rather poor so don't be too upset if you see an unusual patterns of indentations.



#include <iostream>
#include <math.h>
#include <cmath>
#include <fstream>

using namespace std;

int main ()
{
int i,j,k,n,m,x,r,f,s,h,q,t, ans, dim, presin, presout;

cout << "Enter number of rows (at least 3): ";
cin >> n;
cout << "Enter number of columns (at least 3): ";
cin >> m;

//Number of equations developed
dim = (n*m)-(2*n);
s=n-3;

double p[(n*m)+1], b[dim+1];
double eqn[dim+1][(n*m)+1];

for (i=1; i<=n*m; i++)
{
//Fill all nodes with 1 in order to isolate the coeficients of the nodes later
p[i]=1;
}

for (i=1; i<=dim+1; i++)
{
//Fill all nodes with 1 in order to isolate the coeficients of the nodes later
b[i]=0;
}

cout << "Enter imposed inlet pressure: ";
cin >> presin;

for (i=1; i <= (n*m)-(m-1); i+=m)
//Imposed pressures will be over-written
{
p[i] = presin;
}

cout << "Enter imposed outlet pressure: ";
cin >> presout;

for (i=m; i <= n*m; i+=m)
{
p[i] = presout;
p[2] = 1; //enforce p2 = 1...
}

for (i=1, j=1; i<=dim-(m-3) || j<=m*(n-1)+1; i+=m-2, j+=m)
//RHS
{
b[i] = -p[j];
}

for (i=1; i<=dim; i++)
//Fill all equation components with 0 initially
{
for (j=1; j<=n*m; j++)
{
eqn[i][j] = 0;
}
}

//Development of equations
for (i=1; i<=m-2 ; i++)
{
//For top boundary
if (p[i]==1) //If p[i] is a imposed pressure, it skips
{
eqn[i][i] = p[i];
eqn[i][i+1] = - 3*p[i+1];
eqn[i][i+2] = p[i+2];
eqn[i][i+m+1] = p[i+m+1];
}
else
{
eqn[i][i+1] = - 3*p[i+1];
eqn[i][i+2] = p[i+2];
eqn[i][i+m+1] = p[i+m+1];
}
}

for (f=m+1,j=1; f<=(m+1)+(m*(n-3)) || j<=n-2; f+=m, j++)
{
for (i=f; i<=f+(m-3) ; i++)
{
//For the middle
if (p[i]==1)
{
eqn[i-2*j][i-(m-1)] = p[i-(m-1)];
eqn[i-2*j][i+1] = - 4*p[i+1];
eqn[i-2*j][i+2] = p[i+2];
eqn[i-2*j][i+m+1] = p[i+m+1];
eqn[i-2*j][i] = p[i];
}
else
{
eqn[i-2*j][i-(m-1)] = p[i-(m-1)];
eqn[i-2*j][i+1] = - 4*p[i+1];
eqn[i-2*j][i+2] = p[i+2];
eqn[i-2*j][i+m+1] = p[i+m+1];
}
}
}

for (i=dim-(m-3),j=1; i<=dim && j<=m-2 ; i++,j++)
{
//For bottom boundary
if (p[(n-1)*m+j]==1)
{
eqn[i][(n-1)*m+j] = p[(n-1)*m+j];
eqn[i][(n-1)*m+j+1-m] = p[(n-1)*m+j+1-m];
eqn[i][(n-1)*m+j+1] = - 3*p[(n-1)*m+j+1];
eqn[i][(n-1)*m+j+2] = p[(n-1)*m+j+2];
}
else
{
eqn[i][(n-1)*m+j+1-m] = p[(n-1)*m+j+1-m];
eqn[i][(n-1)*m+j+1] = - 3*p[(n-1)*m+j+1];
eqn[i][(n-1)*m+j+2] = p[(n-1)*m+j+2];
}
}

double a[dim+1][dim+1];
int g;

for (i=1; i<=dim; i++)
{
//Fill in the matrix
for (j=1; j<=m-2; j++)
{
a[i][j] = eqn[i][j+1];
}

for (g=m-1, r=0, k=3; g<=(m-1)+((2*n)-m)*s, r<=(n-3)*(m-2), k<=2*n-3; g+=m-2, r+=m-2, k+=2) //issue at a limit
{
for (j=g; j<=(2*m-4)+r; j++)
{
a[i][j] = eqn[i][j+k];
}
}

for (j=(2*m-3)+(n-3)*(m-2); j<=dim; j++)
{
a[i][j] = eqn[i][j+(2*n-1)];
}
}

//Writing the matrix to txt
ofstream myfile ("matrix.txt");
if (myfile.is_open())
{
myfile << dim;
myfile << " 1 ";
myfile << n <<" "<< m <<"\n";
for(i=1;i<=dim;i++)
{
for(j=1;j<=dim;j++)
{
myfile <<a[i][j]<<" ";
}
myfile << "\n";
}

for(j=1;j<=dim;j++)
{
myfile <<b[j]<<" ";
}
}
system ("pause");
return 0;
}



Well I hope this is readable and follow-able.

So if I enter the flow field that is bigger than 30 by 30 (which gives me about 500 unknowns), the computer crashes.
I hope you figure out a way to increase this capacity a bit.

Thanks again.

I've tried the 'static float' bit, but an error apeared as it states that the float is not constant - my matrix size is not fixed. The "easy" way of using static on your matrix and vector declarations works with Standard C++. Standard C++ does not support variable length arrays. (You forgot to tell me what compiler and operating system you are using, but I'm guessing GNU on Windows---maybe Dev-C++ or some such thing---Oh, well, no matter).

Next-to-bottom line: If your compiler supports the non-standard language feature and if it works for you, then who cares about Standards and portability, right?

Bottom line: Your compiler apparently uses the stack for variable length arrays, so that explains why the program chokes, but it also means that you can't make them static. Oh, well...

So: The "next easiest" way might be to use the dvector() and dmatrix() functions in nrutil.c

I don't know if you have access to nrutil.c on your home system, but if you do, then here's a way to get around the stack limitations:

I'll show the lines to change. I commented out the old lines to show what is replacing what (and where).

// Near the top of the file, maybe after the other includes:
#include "nrutil.c"
.
.
.
int main()
{
.
.
.
//double p[(n * m) + 1], b[dim + 1];

double *p;
p = dvector(1, n*m);

double *b;
b = dvector(1, dim);

//double eqn[dim + 1][(n * m) + 1];
double **eqn;
eqn = dmatrix(1, dim, 1, n*m);
.
.
.
//double a[dim + 1][dim + 1];
double **a;
a = dmatrix(1, dim, 1, dim);



Now, if you don't have or can't use nrutil.h, then you will have to learn a little about dynamic allocation in C++. I don't mind showing you if you need it, but first let me know if this works (or can work) for you. No point in a needless waste of board bandwidth (and your personal bandwidth, too), right?


By the way, I can't believe that the following does what you want it to:


for (g=m-1, r=0, k=3;
g<=(m-1)+((2*n)-m)*s, r<=(n-3)*(m-2), k<=2*n-3; //<---wrong, I'm pretty sure
g+=m-2, r+=m-2, k+=2) //issue at a limit


I put the different clauses on different lines so that I could point out what I'm sure is an error. Don't you mean to have "&&" (or, maybe "||"---I didn't try to guess your actual intent) between the conditions. The commas don't connect them logically; the only thing that affects program flow is the last condition.

Regards,

Dave

As you guessed, I am using Dev-C++ at the moment, but eventually I need to introduce some graphics thus perhaps moving to Visual C++ is my intention.

Now, I've made the changes as you shown, and this is showing some big improvements.
It doesn't 'crash' anymore; instead it takes some time and eventually ends up writing the matrix.
At home, I was able to run 35 nodes by 35 nodes (1150 x 1150 matrix), and there's a possibility of achieving something even higher - but I just didn't have the patients.
Then I tested this 1150 x 1150 matrix into the part2 of the code - again, which you've helped me out with - and it took about the same amount of time (or maybe a bit more) and ended up solving it. It took maybe a minute.
So at work, there's a good chance of achieving something greater - but since you say that the system performance has nothing to do with this - maybe not. I'll find out tomorrow and let you know.

And as for the error you pointed out - I didn't notice it until you pointed out, but it's rather bizzarre how it didn't cause any trouble. Nonetheless, I've changed the commas to ||'s.

Thanks so much again, and I'm sure I'll be keeping the threads coming.

So until then,


Pete

So at work, there's a good chance of achieving something greater - but since you say that the system performance has nothing to do with this - maybe not.

Now, don't go putting words into my mouth. What I said was, "the amount of memory available to your C++ program is not directly related to the amount of physical memory in your system."

Performance is another issue. The operating system swaps stuff into and out of physical memory as the various processes require. If there is more physical memory, it just may have better performance since it can keep more of your stuff in physical memory at any given time. It's not guaranteed, but it could happen. I mean, if your process is using a few tens of Megabytes, then having one Gigabyte of physical memory may be just as good as having two Gigabytes. (Or, maybe, not...)

All other things being equal (which processor, processor speed, which Operating System, which compiler, what compiler optimization switches, etc, etc...), more memory is (usually) better.

The three priorities for satisfactory program development:

Priority 1: Don't crash.
Priority 2: Finish before the next ice age.
Priority 3: Get meaningful results


Maybe, just maybe, we have a chance of checking off the first two. Number 3 is the biggie, right? Now to the "Good Stuff."

Regards,

Dave

Footnote: If you move to Visual C++, you will find that variable length arrays are not supported, not even little bittie ones. So maybe code portability is important after all, and getting into the habit of using compiler extensions is, maybe, not a Good Thing. Oh, and by the way, if you are using a really old Visual C++ (version 6, for example), you may find a few tweaks to your standard C++ may be required. No big deal, but something to look forward to, huh?

Hey Dave,

I hope you're still available for some Q&A sessions.
It's that time of my project to start introducing some graphics using openGL, thus I have to switch the platform from Dev C++ to Visual C++ (Visual Basics 2005.
I've just basically copied and pasted the code to the Visual, but the problem seems to occur within the nrutil.c file.
I changed the nrutil.c to nrutil.cpp, but still no use.
I get a bunch of error LNK2005 saying a whole bunch of variables are already defined in the main.

Another source code of mine (part2.exe if you remember), has a similar error when I try to introduce the String Tokenizer function (I bet you know this one). The linker error also developes here.

If you could help me out here, that'd be appricated.

Here's my main:

//HEADER FILES
#include "nrutil.cpp"
#include "matcalc.h" //includes matrix calculation function

using namespace std;

int main ()
{
matcalc();
return 0;
}

And here's my function:

#include <cmath>
#include <fstream>
#include <string>
//#include "nrutil.cpp" I commented this out

//Fill program
//by Peter Kim

using namespace std;

void matcalc ()
{
.
.
.
.
.
code
.
.
.
.
.
cout << "Matrix written to a txt file."<<"\n";
system ("pause");
}

I've just basically copied and pasted the code to the Visual, but the problem seems to occur within the nrutil.c file.

First of all, I just about never actually include C or C++ files in other C or C++ files. I did it in my previous example just to try to get you started. I'll change it a little to be more suitable for projects with multiple files.

Secondly, instead of renaming C files to .cpp in order to get them to play nice together, I just use a command-line switch to do the deed.

Thirdly, I usually use command-line compilation for my command-line programs. Normally I use a Makefile. With Microsoft compilers, there is a utility called nmake that can be used just to recompile the things that change.

In the bin directory of your visual studio installation there should be a file named vcvars32.bat, which sets up environment variables for easiest command-line compiling.

As a matter of fact, there should be a directory called something like "Visual C++ 2005 Express Edition" or some such thing in your Programs stuff that was set up by the installation. Click on "Visual Studio Tools->Visual Studio 2005 Command Prompt"

You can drag that icon to your desktop or start menu, and right-click on it and click "Properties". Under "Start in", enter the path name to where your source files are.

Now, I'll go back to my example. I called it z.cpp, just because I like short names.
Here's the first part of the slightly modified example without the preambulatory comments:


// This is z.cpp
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cstdlib>
#include <cmath>

#include "nr.h"
#include "nrutil.h"
using namespace std;

int main(int argc, char **argv)
{
int i, j, k, m, n, *indx;
.
.
.


Now, suppose I want to compile this along with nrutil.c, ludcmp.c, lubksb.c.

I could do it on a command line. (Kind of messy)
I could do it in a .bat file.
I could do it in a Makefile for use with nmake.

Here's a simple bat file that compiles everything every time. It tells the compiler where to look for include files (nr.h, nrutil.h)
It compiles each source file into a separate obj file and then combines them. After the first time, if you don't change nrutil.c or lubksb.c or ludcmp.c (or the headers) you can comment out those compiler lines in the bat file:


@rem: A little batch file for combining a cpp file
@rem: with Numerical Recipes Version 2 C files
@
@rem: For Microsoft Visual Studio 2005 or 2008
@
@rem: Assumes source file in current directory is z.cpp
@rem: Assumes nrutil.c, nrutil.h and nr.h are in ..\other
@rem: Assumes lubksb.c and ludcmp.c are in ..\recipes
@
@rem: Note you can't combine object files for which some
@rem: were compiled as C and some as C++
@rem:
@rem: Use the command line switch /TP to tell the compiler
@rem: to compile a C file as C++
@rem:
@rem: I made sure that everything is compiled with NRANSI defined
@rem:
@rem: davekw7x
@rem: May, 2008
@
cl -c -I..\other -DNRANSI z.cpp /EHsc
cl -c -I..\other -DNRANSI /TP ..\other\nrutil.c
cl -c -I..\other -DNRANSI /TP ..\recipes\lubksb.c
cl -c -I..\other -DNRANSI /TP ..\recipes\ludcmp.c
cl z.obj nrutil.obj lubksb.obj ludcmp.obj

If you call this m.bat, just enter "m" on the command line. (I like short command lines.)
If you have other files just add lines to compile them.

If they are C++ files, don't forget the /EHsc
If they are C files, don't forget the /TP

If you want to copy nrutil.c, nrutil.h, nr.h, etc. to your working directory, you can leave out the "-I " part and you won't need the "..\..." path stuff. I don't like to do this because it's just too easy to screw up when you have multiple copies of source files sprinkled around the landscape. I just leave everything in its original place.


A much more elegant and practical way for projects with multiple files is to create a Makefile. The nmake program can automatically determine which source files have changed and re-compile only ones that have been updated since the last time. Try the bat file for now and let us know how it works out.


Regards,

Dave

Another source code of mine (part2.exe if you remember), has a similar error

Here's the thing: The Integrated Development Environment is a very useful tool for many things. For simple command-line applications, I use command-line compilation. Why?

1. It's easy to get someone else to compile exactly the same thing in the same way. With Integrated Development Environments there are many, many (many) project settings that may be different from my setup, and it's pretty hard to get them exactly the same without being in the same room and looking over each other's shoulder.

2. If I decide to change compilers, I can see "at a glance" what needs to be compiled and what other stuff needs to be defined instead of somehow magically extracting that information from the depths of project options and who-knows-what from the IDE stuff.

3. I'm funny that way.

So, suppose I have two files: z.cpp and matcalc.cpp. I just compile with


cl z.cpp matcalc.cpp /EHsc


Here's z.cpp:

#include <iostream>

#include "matcalc.h"
using namespace std;

int main ()
{
cout << "***Calling matcalc()***" << endl;
matcalc();
cout << "***Back from matcalc()***" << endl;
return 0;
}

Here's matcalc.h

#include <iostream>
//
// Whatever other includes you need
//
using namespace std;

void matcalc()
{
//
// whatever
//
cout << "Matrix written to a txt file."<<"\n";
cout << "Press 'Enter' to continue . . . ";
cin.get();
}


And matcalc.h can be as simple as

void matcalc();


This gets me exactly what I want. No more, no less. If there are more than two files (or if there are special paths that I need for include or library files) I would probably create a Makefile and use nmake.

If it's a Windows program, it's probably easier with the IDE, since they already have set up things for all of those #include files, #defines, and various library or dll files that are part of a Windows project.

Other times I use the command line, (possibly with a Makefile, or just barefoot as in this example).

And you don't actually need that silly "system(pause)" or equivalent at the end of the program either, since it just returns to the command prompt. Those things (like having to press a key to exit the program) just annoy the hell out of me. I'm funny that way.

If you need nrutil.c (and I guess that you do, for your program), then create a batch file as I showed previously, or a Makefile. The Microsoft nmake program is similar to GNU make (but not as much fun, in my opinion, and not nearly as powerful).

A simple but easily expandable Makefile for nmake might look something like

CXX = cl.exe
LD = cl.exe
DELETE = del

CXXFLAGS= /EHsc
DEFINES = -DNRANSI
INCLUDES= -I..\other

OBJECTS = z.obj nrutil.obj matcalc.obj

TARGET = z.exe

TARGET: $(OBJECTS)
$(LD) z.obj nrutil.obj matcalc.obj

z.obj: z.cpp
$(CXX) -c $(DEFINES) $(INCLUDES) z.cpp $(CXXFLAGS)

nrutil.obj: ..\other\nrutil.c
$(CXX) -c $(DEFINES) $(INCLUDES) ..\other\nrutil.c $(CXXFLAGS)

matcalc.obj: matcalc.cpp
$(CXX) -c $(DEFINES) $(INCLUDES) matcalc.cpp $(CXXFLAGS)

clean:
$(DELETE) *.obj *.exe



Regards,

Dave