Passing zgesv_ a pointer to doubles rather than an array

Open discussion regarding features, bugs, issues, vendors, etc.

Passing zgesv_ a pointer to doubles rather than an array

Postby maberib » Fri Jun 17, 2011 12:56 am

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},
maberib
 
Posts: 3
Joined: Wed May 25, 2011 7:27 pm

Re: Passing zgesv_ a pointer to doubles rather than an array

Postby admin » Wed Jun 22, 2011 8:05 am

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
admin
Site Admin
 
Posts: 490
Joined: Wed Dec 08, 2004 7:07 pm

Re: Passing zgesv_ a pointer to doubles rather than an array

Postby cottrell » Sun Jul 17, 2011 10:01 pm

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
cottrell
 
Posts: 68
Joined: Thu Jan 15, 2009 1:40 pm


Return to User Discussion

Who is online

Users browsing this forum: Bing [Bot] and 2 guests