C-双精度指针与双精度数组之间的细微差别将一个转换为另一个？

分类：编程问答

如何解决C-双精度指针与双精度数组之间的细微差别将一个转换为另一个？

| 我一直在为物理研究项目开发程序。该程序是用C语言编写的，但是使用了fortran函数（它称为\“ zgesv \”，它来自LAPACK和BLAS库）。 这个想法是要解决一个方程组。 LHS.X =向量“ X”的RHS。 int INFO传递给zgesv。如果这些方程无法求解（即LHS是奇异的），则info应该返回值= 0；否则返回0。 通过将double *传递给Solve函数（以下代码中的解决方案1），尝试以“正常”的方式运行程序，即使LHS是单数，INFO也返回0。不仅如此，而且如果我打印出解决方案，那将是数字的灾难-有些小，有些零，有些大。 如果我通过手动创建LHS和RHS     LHS [] = {值1，值2，...};     RHS [] = {值1，值2，...}; 然后将INFO作为预期值返回3，并且解决方案等于RHS（这也是我所期望的）。 如果我声明数组     LHS2 [] = {值1，值2，...};     RHS2 [] = {值1，值2，...}; 并将它们逐个复制到LHS和RHS中，然后将INFO返回为8（对我来说，这与以前的情况不同是很奇怪的。），解决方案等于RHS。 我觉得这肯定是声明数组的两种方式之间的根本区别。我真的无法使用zgesv函数来使muck成为我想要的类型，因为A）它是科学界的标准，并且B）它是用fortran编写的-我还没有从来没有学过。 谁能解释这是怎么回事？是否有一个简单的（最好是计算便宜的）修复程序？我是否应该将我的double *复制到array []中？ 这是我的程序（为测试而修改）： 包括 
#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,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,-1.000000,-0.000000,-3.000000,0.00000};
  double RHS2[] ={0.000000,-2.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\",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\",RHS[i]);
  printf(\"RHS2[%d] = %.16lf\\n\",RHS2[i]);
  printf(\"The difference is %.16lf\",RHS[i] - RHS2[i]);

    }
    i = i + 1;
  }

  printf(\"\\n\");


  double *solution;

  solution = solveUsingLUD(N,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,RHS);

    printf(\"Solution3 = \\n{\");
    for (i = 0; i < 2 * N - 1;)
    {
      printf(\"{%lf + \",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,-1.000000 + 0.000000I,0.000000 + 0.000000I},{0.000000 + 0.000000I,-0.000000 + 0.000000I,-1.000000 + 0.000000I},-0.000000 + -1.000000I,{-1.000000 + 0.000000I,0.000000 + -3.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.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 + -2.000000I},ERROR: info = 8
Solution3 = 
{{0.000000 + 0.000000I},
     

解决方法

        
PI
不能精确地表示为
double
。因此ѭ6都不是。 因此，
sin(kal)
和
cos(kal)
不必分别等于1或0。 （如果要通过将它们分配给\“ int \”来解决此问题，请记住，进行这种转换时C通常会四舍五入。） [编辑] 其实我有一个更好的建议... 使用GCC编译Fortran代码时，您可能要使用
-ffloat-store
。 Fortran代码通常在编写时非常注意双精度数字的量化...但是默认情况下，x86上的中间值可以具有更高的精度（因为浮点单元内部使用80位）。通常，额外的精度只会有所帮助，但是对于某些数字敏感的代码（例如sin（PI / 2）），它可能会导致意外结果。 
-ffloat-store
避免了这种额外的精度。 [编辑2，回应评论] 为了获得更好的正弦和余弦精度，我建议以下几点。 首先，在
main
和其他函数的参数列表中，将
kal
声明为
long double
。 其次，调用
sinl()
和
cosl()
而不是
sin()
和
cos()
。 第三，像这样定义PI： 
#define PI 3.1415926535897932384626433832795029L
将其他所有内容（尤其是
S
和
C
）保留为
double
。看看是否有帮助...     ,        您的代码中有一些语法错误（括号和引号不匹配），但是我想当您在此处复制代码时，这些错误会漏掉。 我的猜测是以下可能会导致问题： 
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： 还有一件事：您正在读取函数getA（）中的
argv[0]
，该函数始终是可执行文件的路径。您应该从
argv[1]
开始阅读。为了安全起见，应检查用户是否提供了足够的参数（
argc - 1 >= N * N
）。     ,首先，我来自Fortran方面，并且没有C经验，因此我所说的内容可能有意义，也可能无关紧要。 BLAS / LAPACK例程希望接收Fortran数组-本质上是连续内存块的第一个元素的地址。假定阵列布局以列为主，并且大小由LDA控制。基本上不会进行任何检查，因此，如果您发送的内容不符合这些期望，那么一切都会变得一团糟。 我从经验中发现（或更确切地说，看到某人这样做并复制了它），从C ++调用LAPACK例程的工作原理如下：如果您使用
std::vector<double>
存储矩阵，并使用
&A[0]
调用LAPACK调用（其中A是
std::vector<double>
），它可以工作。我的理解（也许太天真了）是标准库向量是“有点像” C数组，因此，我想可以将其直接转换为C。 另外，您使用的是复杂的例程，这些例程可能以微妙的方式更改，也可能不会更改。我猜想Fortran的COMPLEX * 16等于
struct{double real_part; double imag_part;}
 最后，可以从netlib处获得LAPACK例程的参考实现，请参见http://www.netlib.org/lapack/complex16/zgesv.f。     ,        这是另一个猜测： 我想知道是否在
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\");
    

版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。

相关推荐

matplotlib报错：AttributeError: module 'backend_interagg' has no attribute 'FigureCanvas'. Did you mean: 'FigureCanvasAgg'?
使用本地python环境可以成功执行 import pandas as pd import matplotlib.pyplot as plt # 设置字体 plt.rcParams[&#39;font.sans-serif&#39;] = [&#39;SimHei&#39;] # 能正确显示负号 p
gitlab登录失败，报错：This challenge page was accidentally cached by an intermediary and is no longer available.
设置时间 控制面板
后端开发常见错误
错误1：Request method ‘DELETE‘ not supported 错误还原：controller层有一个接口，访问该接口时报错：Request method ‘DELETE‘ not supported 错误原因：没有接收到前端传入的参数，修改为如下 参考 错误2：cannot r
docker常见错误
错误1：启动docker镜像时报错：Error response from daemon: driver failed programming external connectivity on endpoint quirky_allen 解决方法：重启docker -&gt; systemctl r
idea常见错误
错误1：private field ‘xxx‘ is never assigned 按Altʾnter快捷键，选择第2项 参考：https://blog.csdn.net/shi_hong_fei_hei/article/details/88814070 错误2：启动时报错，不能找到主启动类 #
pip安装依赖失败
报错如下，通过源不能下载，最后警告pip需升级版本 Requirement already satisfied: pip in c:\users\ychen\appdata\local\programs\python\python310\lib\site-packages (22.0.4) Coll
maven常见错误
错误1：maven打包报错 错误还原：使用maven打包项目时报错如下 [ERROR] Failed to execute goal org.apache.maven.plugins:maven-resources-plugin:3.2.0:resources (default-resources)
cloud项目中常见错误
错误1：服务调用时报错 服务消费者模块assess通过openFeign调用服务提供者模块hires 如下为服务提供者模块hires的控制层接口 @RestController @RequestMapping(&quot;/hires&quot;) public class FeignControl
springboot常见错误
错误1：运行项目后报如下错误 解决方案 报错2：Failed to execute goal org.apache.maven.plugins:maven-compiler-plugin:3.8.1:compile (default-compile) on project sb 解决方案：在pom.
springmvc拦截器中使用redisTemplate报空指针异常
参考 错误原因 过滤器或拦截器在生效时，redisTemplate还没有注入 解决方案：在注入容器时就生效 @Component //项目运行时就注入Spring容器 public class RedisBean { @Resource private RedisTemplate&lt;String
vite报错：npm ERR! command C:WINDOWSsystem32cmd.exe /d /s /c C:UsersychenAppDataLocalTempnpx-dde1c848.cmd
使用vite构建项目报错 C:\Users\ychen\work&gt;npm init @vitejs/app @vitejs/create-app is deprecated, use npm init vite instead C:\Users\ychen\AppData\Local\npm-
VM VirtualBox中安装centos时报错：设置基础软件仓库时出错
参考1 参考2 解决方案 # 点击安装源 协议选择 http:// 路径填写 mirrors.aliyun.com/centos/8.3.2011/BaseOS/x86_64/os URL类型 软件库URL 其他路径 # 版本 7 mirrors.aliyun.com/centos/7/os/x86
kafka常见问题
报错1 [root@slave1 data_mocker]# kafka-console-consumer.sh --bootstrap-server slave1:9092 --topic topic_db [2023-12-19 18:31:12,770] WARN [Consumer clie
hive常见问题
错误1 # 重写数据 hive (edu)&gt; insert overwrite table dwd_trade_cart_add_inc &gt; select data.id, &gt; data.user_id, &gt; data.course_id, &gt; date_format(
装载数据时报错：Failed to execute spark task, with exception 'org.apache.hadoop.hive.ql.metadata.HiveException(Failed to create spark client.)'
错误1 hive (edu)&gt; insert into huanhuan values(1,&#39;haoge&#39;); Query ID = root_20240110071417_fe1517ad-3607-41f4-bdcf-d00b98ac443e Total jobs = 1
flume常见问题
报错1：执行到如下就不执行了，没有显示Successfully registered new MBean. [root@slave1 bin]# /usr/local/software/flume-1.9.0/bin/flume-ng agent -n a1 -c /usr/local/softwa
hadoop集群常见错误
虚拟及没有启动任何服务器查看jps会显示jps，如果没有显示任何东西 [root@slave2 ~]# jps 9647 Jps 解决方案 # 进入/tmp查看 [root@slave1 dfs]# cd /tmp [root@slave1 tmp]# ll 总用量 48 drwxr-xr-x. 2
hive常见错误
报错1 hive&gt; show databases; OK Failed with exception java.io.IOException:java.lang.RuntimeException: Error in configuring object Time taken: 0.474 se
Linux常见错误
报错1 [root@localhost ~]# vim -bash: vim: 未找到命令 安装vim yum -y install vim* # 查看是否安装成功 [root@hadoop01 hadoop]# rpm -qa |grep vim vim-X11-7.4.629-8.el7_9.x
使用sqoop一直卡在：mapreduce.Job: Running job: job_1703173956074_0001
修改hadoop配置 vi /usr/local/software/hadoop-2.9.2/etc/hadoop/yarn-site.xml # 添加如下 &lt;configuration&gt; &lt;property&gt; &lt;name&gt;yarn.nodemanager.res
pythonJavaScriptjavaHTMLPHPreactjsC#AndroidCSSNode.jssqlrpython-3.xMysqLjQueryc++pandasFlutterangularIOSdjangolinuxswifttypescript路由器JSON路由器设置无线路由器h3c华三华三路由器设置华三路由器电脑软件教程arraysdocker软件图文教程Cvue.jslaravelspring-boot