Численное устойчивое прямое замещение

Question

Я должен сделать программу для класса, где программа решает матричную систему с верхней треугольной матрицей.Метод должен быть сделан в C и использовать прямую замену.Затем функция проверяется на некоторые случаи, которые неизвестны, и вы получаете три неудачных / успешных результата.Возвращаемое значение равно 0 для успеха и 1 для числовой ошибки.Однако я не могу заставить программу возвращать 1, когда они этого ожидают.Я пробовал практически на каждом этапе проверять, используются ли nans или infs, и возвращает 1, если они используются.Кроме того, я сделал цикл, чтобы увидеть, дает ли результат исходную правую часть, сделав копию значений, поскольку данная правая часть перезаписана.Я проверяю результат по заданной правой части, вычисляя относительную ошибку.

Тогда у меня вопрос.У кого-нибудь из вас есть идеи относительно того, как я могу провести другие тесты и / или системы, которые бы дали числовую нестабильность?

Если b - правая часть, альфа - это заданная константа, добавленная к диагонали, а R - этоматрица коэффициентов.Там, где элементы в диагонали R плюс альфа не равны нулю.

РЕДАКТИРОВАТЬ: Я сделал следующий код:

#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;

while, в кавычках https://en.cppreference.com/w/c/numeric/math/isnan)

значений NaN никогдасравнивать равные себе или другим значениям NaN.Копирование NaN может изменить его битовую комбинацию.Другой способ проверить, является ли значение с плавающей точкой NaN, состоит в том, чтобы сравнить его с самим собой:

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

Предыдущий тест может быть переписан с использованием библиотечной функции isfinite():

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