Numerical stable forward substitution

Question

I have to make a program for class where the program solves a matrix system with an upper triangular matrix. The method has to be made in C and use forward substitution. The function is then tested against some cases which are not known and you get three failed/succeeded results back. The return value is 0 for success and 1 for numerical failure. However, I cannot get the program to return 1 when they expect it. I have tried testing, at pretty much every step, if nans or infs are in use and return 1 if they are. Furthermore, I have made a loop to see if the result gives the original righthand side, by making a copy of the values as the given right-hand side is overwritten. I test the result against the given right-hand side by calculating the relative error.

My question is then. Do any of you have ideas as to how I can make other tests and/or systems, which would give numerical instability?

If b is the righthand side, alpha is a given constant added to the diagonal and R is the coefficient matrix. Where the elements in the diagonal of R plus alpha is not zero.

EDIT: The code I made is the following

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

And the function

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

Since you're not really asking any C language related questions, why not drop the C tag and add Algorithm? — jwdonahue, Sep 28 '18 at 18:37
Perhaps you have a programming error rather than needing to devise more tests. Make sure you have full compiler warnings enabled. — Weather Vane, Sep 28 '18 at 18:46
Please note that conditions like `err == NAN` are always false. See e.g. https://stackoverflow.com/a/1923933/4944425 . Your function is also returning -1, not 1 (as stated in the question) in case of failure, by the way. — Bob__, Sep 28 '18 at 21:13

score 0 · Accepted Answer · answered Sep 28 '18 at 21:38

One error in OP's code is directly comparing a double value with NAN, like in

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

while, quoting 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; }

The previous test could be rewritten using the library function isfinite():

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

Numerical stable forward substitution

1 Answers1