C - незначительные различия между указателями на двойники и массивами двойников. Преобразовать одно в другое? - PullRequest
Новая
       11

C - незначительные различия между указателями на двойники и массивами двойников. Преобразовать одно в другое?

1 голос
/ 18 июня 2011

Я работал над программой для исследовательского проекта по физике. Программа написана на C, но использует функцию fortran (она называется "zgesv" и из библиотек LAPACK и BLAS).

Идея состоит в том, чтобы решить систему уравнений. LHS.X = RHS для вектора "X". int INFO передается в zgesv. Если уравнения не могут быть решены (то есть LHS является единственным), информация должна возвращать значение! = 0;

Пытаясь запустить мою программу как "нормальную", передав мой double * в функцию решения (решение 1, в следующем коде), INFO возвращается как 0 - даже если LHS единственное число. И не только это, но если я распечатаю решение , это катастрофа чисел - некоторые маленькие, некоторые ноль, некоторые огромные.

Если я создаю LHS и RHS вручную, выполнив LHS [] = {значение 1, значение 2, ...}; RHS [] = {значение 1, значение 2, ...}; Затем INFO возвращается как 3, как ожидалось, и решение равно RHS (что также является тем, что я ожидал).

Если я объявлю массивы LHS2 [] = {значение 1, значение 2, ...}; RHS2 [] = {значение 1, значение 2, ...}; и копировать их в LHS и RHS элемент за элементом, затем INFO возвращается как 8 (мне странно, что он отличается от предыдущего случая.), и решение равно RHS .

Мне кажется, это должно быть фундаментальное различие между двумя способами объявления массива. На самом деле у меня нет доступа к гадости с функцией zgesv, чтобы она принимала типы, которые я хочу, так как A) это стандарт в научном сообществе, и B) он написан на фортране - который я никогда не изучал.

Может кто-нибудь объяснить, что здесь происходит? Есть ли простое (и желательно дешевое) исправление? Должен ли я просто скопировать свой double * в массив []?

Вот моя программа (модифицированная для тестирования):

включает 

#include <stdlib.h>
#include <math.h>
#define PI 3.1415926535897932384626433832795029L

#define ERROR_VALUE 911.911


int* getA(int N, char* argv[])
{
  int i;
  int* AMatrix;
  AMatrix = malloc(N * N * sizeof(int));

  if (AMatrix == NULL)
  {
    printf("Failed to allocate memory for AMatrix. Exiting.");
    exit (EXIT_FAILURE);
  }

  for (i = 0; i < N * N; i++)
  {
    AMatrix[i] = atoi(argv[i + 1]);
  }

  return AMatrix;

}

double* generateLHS(int N, int* AMatrix, int TAPs[], long double kal)
{

  double S, C;
  S = sinl(kal);
  C = cosl(kal);
  printf("According to C, Sin(Pi/2) = %.25lf and Cos(Pi/2) = %.25lf", S, C);
  //       S = 1;
  //       C = 0;
  double* LHS;
  LHS = malloc(N * N * 2 * sizeof(double));

  if (LHS == NULL)
  {
    printf("Failed to allocate memory for LHS. Exiting.");
    exit (EXIT_FAILURE);
  }

  int i;
  for (i = 0; i < N * N; i++)
  {
    LHS[2 * i] = -1 * AMatrix[i];  
    LHS[(2 * i) + 1] = 0;
  }


  for (i = 0; i <= 2 * N * N - 2; i = i + (2 * N) + 2)
  {
    LHS[i] = LHS[i] + (2 * C); 
  }

  int j;
  for (i = 0; i <= 3; i++)
  {
    j = 2 * N * TAPs[i] + 2 * TAPs[i];
    LHS[j] = LHS[j] - C;
    LHS[j + 1] = LHS[j + 1] - S;
  }

  return LHS;
}

double* generateRHS(int N, int inputTailAttachmentPoint, long double kal)
{

  double* RHS;
  RHS = malloc(2 * N * sizeof(double));

  int i;
  for (i = 0; i < 2 * N; i++)
  {
    RHS[i] = 0.0;
  }
  RHS[2 * inputTailAttachmentPoint + 1] = - 2 * sin(kal);
  return RHS; 
}


double* solveUsingLUD(int N, double* LHS, double* RHS)
{
  int INFO; /*Info is changed by ZGELSD to 0 if the computation was carried out successfully. Else it changes to some none-zero integer. */
  int ione = 1;
  int LDA = N;
  int LDB = N;
  int n = N;
  int* IPV = malloc(N * sizeof(int));

  if (IPV == NULL)
  {
    printf("Failed to allocate memory for IPV. Exiting.");
    exit (EXIT_FAILURE);
  }

  zgesv_(&n, &ione, LHS, &LDA, IPV, RHS, &LDB, &INFO);
  free(IPV);

  if (INFO != 0)
  {

    printf("\n ERROR: info = %d\n", INFO); 
  }
  return RHS;
}


void printComplexVectors(int numberOfRows, double* matrix)
{
  int i;

  for (i = 0; i < 2 * numberOfRows - 1; i = i + 2)
  {
    printf("%f + %f*i    \n", matrix[i], matrix[i + 1]);
  }
  printf("\n");
}


int main(int argc, char* argv[])
{
  int N = 8;
  int* AMatrix;
  AMatrix = getA(N, argv);
  int TAPs[]={4,4,4,3};
  long double kal = PI/2;




  double *LHS, *RHS;
  LHS = generateLHS(N, AMatrix, TAPs,kal);
  int i;
  RHS = generateRHS(N, TAPs[0],kal);





  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("}");


  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,0.00000};
  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};

  printf("comparing LHS and LHS2\n");
  for (i = 0; i < 2 * N * N;)
  {
    if (LHS[i] != LHS2[i]) {
      printf( "LHS difference at index %d\n", i);
      printf("LHS[%d] = %.16lf\n", i, LHS[i]);
  printf("LHS2[%d] = %.16lf\n", i, LHS2[i]);
  printf("The difference is %.16lf\n", LHS[i] - LHS2[i]);
    }
    i = i + 1;
  }

  printf("\n");
  printf("comparing RHS and RHS2\n");
  for (i = 0; i < 2 * N;)
  {
    if (RHS[i] != RHS2[i]) {
      printf( "RHS difference at index %d\n", i);
      printf("RHS[%d] = %.16lf\n", i, RHS[i]);
  printf("RHS2[%d] = %.16lf\n", i, RHS2[i]);
  printf("The difference is %.16lf", RHS[i] - RHS2[i]);

    }
    i = i + 1;
  }

  printf("\n");


  double *solution;

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


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


    for (i = 0; i < 2 * N * N;)
    {
      LHS[i] = LHS2[i];
      i = i + 1;
    }


    for (i = 0; i < 2 * N;)
    {
      RHS[i] = RHS2[i];
      i = i + 1;
    }
    solution = solveUsingLUD(N,LHS,RHS);

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

Я использую строку компиляции 

gcc -lm -llapack -lblas PiecesOfCprogarm.c -Wall -g

и я выполняю с: 

./a.out 0 0 0 0 1 1 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 1 1 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0

Что дает на выходе 

According to C, Sin(Pi/2) = 1.0000000000000000000000000 and Cos(Pi/2) = -0.0000000000000000000271051
 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},
{}comparing LHS and LHS2
LHS difference at index 0
LHS[0] = -0.0000000000000000
LHS2[0] = 0.0000000000000000
The difference is -0.0000000000000000
LHS difference at index 18
LHS[18] = -0.0000000000000000
LHS2[18] = 0.0000000000000000
The difference is -0.0000000000000000
LHS difference at index 36
LHS[36] = -0.0000000000000000
LHS2[36] = 0.0000000000000000
The difference is -0.0000000000000000
LHS difference at index 54
LHS[54] = -0.0000000000000000
LHS2[54] = 0.0000000000000000
The difference is -0.0000000000000000
LHS difference at index 72
LHS[72] = 0.0000000000000000
LHS2[72] = -0.0000000000000000
The difference is 0.0000000000000000
LHS difference at index 90
LHS[90] = -0.0000000000000000
LHS2[90] = 0.0000000000000000
The difference is -0.0000000000000000
LHS difference at index 108
LHS[108] = -0.0000000000000000
LHS2[108] = 0.0000000000000000
The difference is -0.0000000000000000
LHS difference at index 126
LHS[126] = -0.0000000000000000
LHS2[126] = 0.0000000000000000
The difference is -0.0000000000000000

comparing RHS and RHS2


 Solution = 
{{1.0000000000000000 + -0.0000000000000000I},
{-1.0000000000000000 + -0.0000000000000000I},
{-0.0000000000000000 + 0.0000000000000000I},
{-0.0000000000000000 + 0.0000000000000000I},
{0.6000000000000000 + 0.2000000000000000I},
{-0.0000000000000000 + -0.0000000000000000I},
{-6854258945071195.0000000000000000 + 4042255275298396.0000000000000000I},
{6854258945071195.0000000000000000 + -4042255275298396.0000000000000000I},

 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},

 ERROR: info = 8
Solution3 = 
{{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},

Ответы [ 4 ]

3 голосов
/ 18 июня 2011

У вас есть некоторые синтаксические ошибки в вашем коде (несоответствующие скобки и кавычки), но я предполагаю, что они возникли, когда вы скопировали код здесь.

Я предполагаю, что следующее может вызвать проблему: 

solveSystemByLUD(int N, double *LHS, double* RHS, ...)
{
...
}

Это объявляет функцию, которая возвращает int, но вы возвращаете double *.Добавьте тип возврата и посмотрите, что изменится.

РЕДАКТИРОВАТЬ: После вашего разъяснения все еще есть некоторые недостатки - код по-прежнему не может скомпилироваться из-за этой строки: 

   * SUBROUTINE ZGESV( N, Nrhs, A, LDA, IPIV, B, LDB, INFO ) Solves A * X = B and stores the answer in B.   If the solver worked, INFO is set to 0.
   */

Отсутствует вводный комментарий.

Также вы присваиваете результаты sin ()и cos () - целочисленные переменные (по умолчанию C и S - целые числа): 

  S = sin(kal);
  C = cos(kal);

Судя по вашему коду, это не то, что вам нужно.

Последнее, на что нужно обратить внимание, это то, чтов LHS2 отсутствует последний элемент (sizeof(LHS2)/sizeof(double) равно 127 вместо 2 * 8 * 8 = 128).Это означает, что вы читаете и пишете за пределами массива, что вызывает неопределенное поведение и может вызвать проблемы, которые вы видите.

РЕДАКТИРОВАТЬ 2: Еще одна вещь: вы читаете argv[0] в функции getA (), которая всегда является путем к исполняемому файлу.Вы должны начать читать с argv[1].И чтобы быть в безопасности, вы должны проверить, предоставил ли пользователь достаточно аргументов (argc - 1 >= N * N).

1 голос
/ 19 июня 2011

Прежде всего: я из Фортрана и у меня нет опыта C, поэтому то, что я говорю, может иметь или не иметь значение.

Подпрограммы BLAS / LAPACK ожидают получения массивов Фортрана - по существу, это адрес первого элемента непрерывного фрагмента памяти.Предполагается, что макет массива должен быть главным столбцом, а размер контролируется LDA.По сути, никаких проверок не производится, поэтому, если то, что вы отправляете, не соответствует этим ожиданиям, весь ад проваливается.

То, что я обнаружил из опыта (или, скорее, увидел, что кто-то делает это и скопировал его) при вызове подпрограмм LAPACK из C ++, заключается в следующем: если вы используете std::vector<double> для хранения своих матриц,и вызывать вызовы LAPACK, используя &A[0] (где A - std::vector<double>), это работает.Мое (может быть, слишком наивное) понимание состоит в том, что стандартный вектор библиотеки является "чем-то вроде" массива C, так что, я думаю, это может быть напрямую переведено в C.

Кроме того, вы используете сложные подпрограммы, который может или не может изменить что-то тонким образом.Я предполагаю, что КОМПЛЕКС Фортрана * 16 эквивалентен struct{double real_part; double imag_part;}

Наконец, эталонная реализация подпрограмм LAPACK легко доступна из netlib, здесь http://www.netlib.org/lapack/complex16/zgesv.f

0 голосов
/ 19 июня 2011

PI не совсем точно представляется как double.Поэтому ни PI/2.

Поэтому sin(kal) и cos(kal) не обязательно равны 1 или 0 соответственно.

(Если вы пытаетесь это исправить, назначив их"int", помните, что C обычно округляет вниз при выполнении такого преобразования.)

[править]

На самом деле у меня есть лучшее предложение ...

При компиляции кода на Фортране с помощью GCC вы, вероятно, захотите использовать -ffloat-store.Код на Фортране часто пишется с внимательным вниманием к квантованию чисел с двойной точностью ... Но по умолчанию промежуточные значения в x86 могут иметь дополнительную точность (поскольку модуль с плавающей запятой использует 80 бит внутри себя).Обычно дополнительная точность только помогает, но для некоторого кода, чувствительного к числовому значению (например, sin (PI / 2)), это может привести к неожиданным результатам.

-ffloat-store позволяет избежать этой дополнительной точности.

[править 2, в ответ на комментарии]

Чтобы получить лучшую точность от sin и cos, я предлагаю следующее.

Сначала объявим kal как long double,и в main, и в списках аргументов для других функций.

Во-вторых, вызовите sinl() и cosl() вместо sin() и cos().

В-третьих, определитеPI вот так: 

#define PI 3.1415926535897932384626433832795029L

Оставьте все остальное (особенно S и C) как double.Посмотрите, поможет ли это ...

0 голосов
/ 18 июня 2011

Вот еще одно предположение:

Интересно, если внутри generateLHS() и / или generateRHS() вы генерируете double значений с очень маленькими ошибками, которые отображаются как ноль (или, по крайней мере, ноль в дробной части), когда вы отображаете их с printf(), но различия влияют на расчеты в zgesv_().

Можете ли вы попробовать добавить эти циклы в тестовый прогон сразу после объявлений LHS2 и RHS2: 

printf("comparing LHS and LHS2\n");
for (i = 0; i < 2 * N * N;)
{
    if (LHS[i] != LHS2[i]) {
        printf( "LHS difference at index %d\n", i);
    }
    i = i + 1;
}
printf("\n");
printf("comparing RHS and RHS2\n");
for (i = 0; i < 2 * N;)
{
    if (RHS[i] != RHS2[i]) {
        printf( "RHS difference at index %d\n", i);
    }
    i = i + 1;
}
printf("\n");
...