数值稳定正向替换

Question

我必须为 class 编写一个程序，该程序求解具有上三角矩阵的矩阵系统。该方法必须用 C 语言编写并使用正向替换。然后针对一些未知的情况对该函数进行测试，您会得到三个 failed/succeeded 结果。 return 值为 0 表示成功，1 表示数值失败。但是，当他们期望时，我无法将程序设置为 return 1 。我几乎在每一步都尝试测试是否正在使用 nans 或 infs，如果正在使用 return 1。此外，我做了一个循环来查看结果是否给出了原始的右侧，方法是在给定的右侧被覆盖时复制值。我通过计算相对误差来测试给定右侧的结果。

我的问题是。你们有没有人知道我如何进行其他测试 and/or 系统，这会导致数值不稳定？

如果 b 是右侧，则 alpha 是添加到对角线上的给定常数，R 是系数矩阵。其中R加上alpha的对角线元素不为零。

编辑： 我做的代码如下

#include <math.h>
#include <stdlib.h>
#include <string.h>

int fwsubst(unsigned long n,double alpha,double **R,double *b){
double *x = malloc(sizeof(double)*n);
double row, err = 0,t,tmp;

memcpy(x,b,sizeof(double)*n);

for (size_t k = 0; k < n; k++) {
tmp = (b[k]-summation(R,b,k))/(R[k][k]+alpha);
  if(tmp ==-NAN || tmp ==NAN || tmp ==INFINITY || tmp ==-INFINITY)
    return -1;
  b[k]=tmp;
} 

for (size_t k = 0; k < n; k++) {
  row = 0;
  for (size_t i = 0; i < k; i++) {
    row += R[i][k] * b[i];
  }
  row += b[k]*(R[k][k]+alpha);
  if (row==NAN || row==INFINITY)
    return -1;

  if(x[k]==0)
    t = fabs(row);
  else
    t = fabs(1-row/x[k]);

  err += t;
}

if(err>1e-15 || err==-NAN || err==NAN || err==INFINITY || err==-INFINITY)
  return -1;

return 0;
}

和函数

double summation(double **R, double* b, size_t k){
double sum=0;
for (size_t i = 0; i < k; i++) {
  sum += R[i][k]*b[i];
}

return sum;
}

Answer 1

OP 代码中的一个错误是直接将 double 值与 NAN 进行比较，如

if(tmp ==-NAN || tmp ==NAN || tmp ==INFINITY || tmp ==-INFINITY)
    return -1;

同时引用 https://en.cppreference.com/w/c/numeric/math/isnan)

NaN values never compare equal to themselves or to other NaN values. Copying a NaN may change its bit pattern. Another way to test if a floating-point value is NaN is to compare it with itself:

bool is_nan(double x) { return x != x; }

可以使用库函数重写之前的测试isfinite():

if ( !isfinite(tmp) )
    return -1;

数值稳定正向替换

Numerical stable forward substitution

c

algorithm

equation

numeric

matrix