Passing an LAPACK (Fortran linear algebra library) function a pointer to double rather than an array of doubles

maberib (2)
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. For those that don't know: LAPACK is a linear algebra library written in Fortran. My code is in C and compiled using gcc.

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: 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
 
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:
1
2
3
4
5
 
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:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
 

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:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
 
 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:
1
2
3
4
5
6
7
8
9
10
11
12
13
 
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:
1
2
3
4
5
6
7
8
9
10
11
 
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},
Last edited on
Topic archived. No new replies allowed.