Passing zgesv_ a pointer to doubles rather than an array

[maberib]

Hey guys,

I'm hoping one of you can help me out with this.

F.Y.I. The problem might be solvable with just the first paragraph and first code block, so you may not have to read this whole thing.

I have a big list of matrix equations that I need to solve (A.X= B, solve for X). I use LAPACK's LU decomposition routine, zgesv, to solve them. Sometimes, A will be singular. One of zgesvs arguments, int INFO, is supposed to return non zero if it's singular.

When I pass it a pointer that was set equal to a malloc'd array, it does not detect the singularity and returns INFO = 0 with a garbage solution. But when I pass it the same values stored in a "normal" array, it returns INFO = 3, as expected. I think there's a subtlety that I'm missing here. Is what I'm doing fundamentally wrong in anyway? I'd rather not copy my malloc'd array element by element into a "normal" array every time.

For example, in pseudocode:

Code: Select all

double *ReturnMatrix1(int N, int* b)
{
double *NewMatrix = malloc(N * N * sizeof(doubles));
 ... //do some stuff to fill in NewMatrix entry by entry
 ...
return NewMatrix;
}

double *ReturnMatrix2(int N, int* b, int c[])
{
double *NewMatrix = malloc(N * N * sizeof(doubles));
 ... //do some stuff to fill in NewMatrix entry by entry
 ...
return NewMatrix;
}

solveSystem(int N, double *LHS, double *RHS)
{
int ione = 1;
int* IPV = malloc(N * sizeof(int));
int INFO;
zgesv_(&N, &ione, LHS, &N, IPV, RHS, &N, &INFO);
return RHS;
}
int main()
{
  int N = 8;
  double *LHS, *RHS;
  int *b;
  int c[] = {1,2,3,4};
...
// do some stuff to fill in b entry by entry
...
LHS = ReturnMatrix1(N, b);
RHS = ReturnMatrix2(N, b, c);
RHS = SolveSystem(N, RHS, LHS);



Now RHS contains garbage. The last two entries of RHS are very large integer numbers.

The matrices that are causing me trouble are:

Code: Select all

LHS:
{0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-0.000000,-3.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000};

RHS:
{0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-2.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000};


I've tested that LHS and RHS are actually holding what I want by using the code:

Code: Select all


double* solveUsingLUD(int N, double* LHS, double* RHS)
{
  int INFO;
  int ione = 1;
 
  int* IPV = malloc(N * sizeof(int));
 
  if (IPV == NULL)
  {
    printf("Failed to allocate memory for IPV. Exiting.");
    exit (EXIT_FAILURE);
  }

  zgesv_(&N, &ione, LHS, &N, IPV, RHS, &N, &INFO);
  free(IPV);
 
  if (INFO != 0)
  {
   
    printf("\n ERROR: info = %d\n", INFO);
  }
  return RHS;
}

int main()
{
...
...
printf("\n LHS = \n{{");
 for (i = 0; i < 2 * N * N - 1;)
 {
   printf("%lf + ", LHS[i]);
   i = i + 1;
   printf("%lfI", LHS[i]);
   i = i + 1;
   
   if ((int)(.5 * i) % N == 0)
   {
      printf("},\n{");
   }
   else
   {
     printf(",");
   }
 }
 printf("}");

 printf("\n RHS = \n{");
  for (i = 0; i < 2 * N - 1;)
 {
   printf("{%lf + ", RHS[i]);
   i = i + 1;
   printf("%lfI}, ", RHS[i]);
   i = i + 1;

    printf("\n");
   
 }
 printf("}");
 RHS = solveUsingLUD(N, LHS, RHS);
 printf("\n Solution = \n{");
  for (i = 0; i < 2 * N - 1;)
 {
   printf("{%lf + ", RHS[i]);
   i = i + 1;
   printf("%lfI},", RHS[i]);
   i = i + 1;
   printf("\n");
 }
}


Which gives output:

Code: Select all

 LHS =
{{0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I},
{0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I},
{0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I},
{0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + -1.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I},
{-1.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,-0.000000 + -3.000000I,0.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I},
{-1.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I},
{0.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I},
{0.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I,-1.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I,0.000000 + 0.000000I},
{}
 RHS =
{{0.000000 + 0.000000I},
{0.000000 + 0.000000I},
{0.000000 + 0.000000I},
{0.000000 + 0.000000I},
{0.000000 + -2.000000I},
{0.000000 + 0.000000I},
{0.000000 + 0.000000I},
{0.000000 + 0.000000I},
}
 Solution =
{{4.000000 + 4.000000I},
{-4.000000 + -4.000000I},
{0.000000 + 0.000000I},
{0.000000 + 4.000000I},
{-0.800000 + 2.400000I},
{-4.000000 + 0.000000I},
{33626877217699712.000000 + 4803839602528530.000000I},
{-33626877217699708.000000 + -4803839602528530.000000I},
}


where as when I run the code:

Code: Select all

double LHS2[] = {0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-0.000000,-3.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,-1.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000};
double RHS2[] ={0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-2.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000};
RHS = solveUsingLUD(N, LHS2, RHS2);

 printf("\n Solution2 = \n{");
  for (i = 0; i < 2 * N - 1;)
 {
   printf("{%lf + ", RHS[i]);
   i = i + 1;
   printf("%lfI},", RHS[i]);
   i = i + 1;
   printf("\n");
 }


I get what I expect:

Code: Select all

ERROR: info = 3

 Solution2 =
{{-0.000000 + 0.000000I},
{0.000000 + 0.000000I},
{0.000000 + 0.000000I},
{0.000000 + 0.000000I},
{0.000000 + -2.000000I},
{0.000000 + 0.000000I},
{0.000000 + 0.000000I},
{0.000000 + 0.000000I},


[admin]

Hi,
I don't see anything wrong in the pieces of code. If you provide me a complete code that isolate the problem, I could take a further look.
One advice would be to use the Standard C language APIs for LAPACK
(http://www.netlib.org/lapack/#_standard ... for_lapack).
It would make your code simpler and more readable.
Julie

[cottrell]

I can attest from experience with many lapack functions that it
makes no difference whether you pass a pointer to correctly
malloc'd memory or the name of an array on the stack; in terms
of the C programming language it would be truly weird if there
were any difference.

I note that there's an error in the snippets you posted: namely
the use of "sizeof(doubles)" rather than "sizeof(double)". I presume
this is not what you actually have in your code, or it would not compile.

The most likely explanation is an error in your code. I recommend
compiling your program with debugging symbols enabled and
running it under valgrind [ http://valgrind.org/ ].

Allin Cottrell