Hi,

I'm using NR to solve a linear system - but the size of the system increases with time. I'm using memory allocation as sugested by davekwx7 back in May.

Basically the matrix size is determined then the calculation is carried, and this is looped so that it repeates the process of scan for size, then calculate.

I get the allocation failure 2 in matrix() error when this matrix becomes too big; 288 by 288 to be exact. I actually need this to be able to handle something a lot bigger than that, perhaps up to a 500 by 500.

Is there any way of improving this?
I have a 3.0 Gb ram, if that is any relevent, and a P4 3.4 GHz.

Here's my code:

//HEADER FILES
#include <iostream>
#include <iomanip>
#include <math.h>
#include <cmath>
#include <fstream>
#include <string>
#include <stdlib.h>
#include <cstdlib>

//Matrix create code
//by Peter Kim
//writes "matrix.txt" and "solution.txt"

using namespace std;

int obj=0;

void matcalc3 ()
{

//Variables
char read;
int i, j, k;
int ctner=0, fin=0;

ifstream inData;
string line;
char * name = "test.txt";

inData.open(name);

if (!inData) {
cerr << "Unable to open file";
exit(1); // call system to stop
}

if (inData.is_open())
{
//Reads first line of the txt (#row, #col, #unknown)
inData >> n;
inData >> m;
inData >> dim;


eqn = dmatrix(1, dim, 1, n*m);

for (i=1; i<=dim; i++) {
for (j=1; j<=n*m; j++) {
eqn[i][j] = 0;
}
}


b = dvector(1, dim);

for (i=1; i<=dim; i++)
b[i] = 0;
.
.
.
.
//some irrelevent stuff here
inData.close();
.
.
.
//more irrelevent stuff here
for (i=1; i<=dim; i++) {
if (b[i] == -0)
b[i] = 0;
}


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


//double **a;
a = dmatrix(1, dim, 1, dim);

for(i=1;i<=dim;i++) {

for(j=1;j<=dim;j++) {

a[i][j] = 0;
}

}

int clm = 1;

for (i=1; i<=dim; i++) { //Fill in the matrix
for (k=1; k<=n*m; k++)
{
if (type[k]==2)
{
a[i][clm] = eqn[i][k];
clm +=1;
}
}
clm = 1;
}



//Writing the matrix to txt
ofstream myfile ("matrix.txt");
if (myfile.is_open())
{
myfile << dim;
myfile << " 1 "; //number of matrix
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]<<" ";
}
}
//cout << "Matrix written."<<"\n";

//} //end of the else if (fin != 3)



//NR stuff

int *indx;
int num, job=0, job2=0;
double p, **aa;
double *x;
double sum;

indx = ivector(1, dim); /* holds pivot info from ludcmp for lubskp*/
aa = matrix(1, dim, 1, dim); /* a copy of the matrix to send to ludcmp */
x = vector(1, dim); /* Copy rhs to x and send to lubksp */

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

//cout << scientific << showpos;

for (i=1; i<=dim; i++) {
for (j=1; j<=dim; j++) {
//cout << "a[" << i <<"]["<<j<<"]="<<a[i][j]<<"\n";
aa[i][j] = a[i][j];
}
}

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

/* Solve equations for each right-hand vector */
for (k = 1; k <= dim; k++) {
x[k] = b[k]; /* save b, work with x */
}

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

//cout << "Solution:" << endl;

for (j=0; j<=n-1; j++)
{

for (i=1; i<=m; i++)
{
if (type[j*m+i]==0) {
//cout << setw(10) << setprecision (3) << pres[j*m+i];
}
else if (type[j*m+i]==3) {
//cout << setw(10) << " ";
}
else if (type[j*m+i]==2) {
job+=1;
//cout << setw(10) << setprecision (3) << x[job];
pres[j*m+i] = x[job];
}
else if (type[j*m+i]==1) {
//cout << setw(10) << setprecision (3) << pres[j*m+i]; //used to be just 'i'
}

else if (type[j*m+i]==4) {
//cout << setw(10) << setprecision (3) << "#"; //used to be just 'i'
}



}
//cout << '\n';

}

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

//Writing solution to txt
ofstream solfile ("solution.txt");
if (solfile.is_open())
{
solfile << "Solution:" <<"\n";
for (j = 0; j <= n-1; j++)
{
for (i=1; i<=m; i++)
{
if (type[j*m+i]==0) {
solfile << setw(10) << setprecision (3) << pres[j*m+i];
}
else if (type[j*m+i]==3) {
solfile << setw(10) << " ";
}
else if (type[j*m+i]==2) {
job2+=1;
solfile << setw(10) << setprecision (3) << x[job2];
}
else if (type[j*m+i]==1) {
solfile << setw(10) << setprecision (3) << pres[j*m+i];
}
else if (type[j*m+i]==4) {
solfile << setw(10) << setprecision (3) << "#"; //used to be just 'i'
}
}
solfile << "\n";
solfile << "\n";
}
}

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

//cout << "Solution file written." << "\n";

//cout << "xcount = " << xcount << "\n";


} //eop


There you have it...
If you need to know anything else, feel free to post replies.
If there's a way to avoid using memory allocation, WITHOUT destroying the structure so much, that would also be welcomed to know.
At this stage of my grad studies, I can not afford to start from scartch.

Thank you!

Keep in mind that above is only a header file from my entire code, I've got a bit more stuff then that.
So you won't be able to complie it... and I know I'm sorry for that.
But I just can't give away too much as this is going into my thesis...

So you won't be able to complie it...That's putting it mildly.
..But I just can't give away too much ...
That's entirely understandable.

Here's the thing: You need someone who can look at your code (the whole enchilada) to look for memory problems. There are tools that may be able to spot leaks (valgrind comes to mind), but they all have their little quirks,and there is a learning curve (primarily false positives) that sometimes takes a little time. (But we all have our little quirks, right?)

Here's my recommendation: Look at each and every block of code where you allocate storage: matrix, dmatrix, vector, ivector, etc.

Make absolutely sure that you deallocate that storage when you are through with it. If the block is inside a function then unless you are returning a pointer to that storage, to the calling program, you should free up the storage before you leave that block. If you do return a pointer to that storage, then make absolutely sure the calling function frees it up when it is through. If the calling function returns that pointer value to some higher layer routine, then the function that called that one will be responsible for freeing the storage.

For example in your code, I see the following:

eqn = dmatrix(1, dim, 1, n*m);


I don't see where you free up that storage and I don't see where you pass the value to some other routine that can free it up when it is through. Maybe that is in the code that you couldn't show us.

Anyhow, that's the idea. Maybe you can follow through but it's quite impossible for us to guess how many ways there are go go wrong. (Well, maybe not impossible for some, but pretty much beyond my personal psychic abilities.)

Regards,

Dave

Hi Dave,

Thanks for the tip once again...
I'm not at work at the moment so I can't really try out your suggestion, but I get the basic idea.
As you spotted, I don't think I free up the space for the for eqn - matter of fact I've got a few more dvector like this which I don't free up the space after.
I only free up the ones you see at the bottom of my sourced code:
b,a,aa,x,indx

If I free up the space for other dmatrix and dvectors, that should help out a bit, right?

So where do I exactly add the free_vector and free_matrix?
At the end of the usage as it is for b,a,aa,x,and indx, right?
So if the program updates the matrix after it solves it, I should free it up right after it solves it, so that the memory is freed BEFORE it creates the matrix... correct?

I can't wait to go back in the lab and try this out!
Thanks Dave always.

Dave,

I tried what you suggested, but it seems to worsen the problem, as either it crashes eariler, or bugs the program.

At this point, I think someone needs to look at my whole enchilada as you mentioned, and I can't really think of anyone but you.

However, I do not want to post it up here so that some smart ass could jack it from me... (it's not a nobel prize worthy work, but I do want to secure it).

So is it possible for me to send to through your email perhaps? Will you be willing to look at my whole code and correct this memory problem? I would not know how to thank you if you do that.

If you are willing, please leave me your email.

Thanks.

Nevermind Dave!

I freed up the memory for eqn only - which is the one used in that particular function.
The rest I left them as they are, and the results are looking better.
It solved up to 520 by 520, and it looks to be capable or doing even more - i just ran out of time to try out a bigger one.

Thanks always!
Pete