c - 雅可比方法适用于 double 型，但适用于浮点型。怎么了？

Question

我编写了一个小程序，用雅可比(迭代)方法求解 n 个方程组。下面是代码:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
int main() {

float *a, *b, *x, *xnew, temp;
int i, j, k, maxiter=10000000, n=4;

a = malloc(n*n*sizeof(*a));
b = malloc(n*sizeof(*b));
x = malloc(n*sizeof(*x));
xnew = malloc(n*sizeof(*xnew));

srand((unsigned) time(NULL));

//  Filling the matrix
for (i=0;i<=n-1;i++) {
    for (j=0;j<=n-1;j++) {
        a[n*i+j] = rand()%60;
    }
    b[i] = rand();
    x[i] = rand();
    xorg[i]=x[i];
}

//  Establishing diagonal dominance
for (i=0;i<=n-1;i++) {
    temp=0;
    for (j=0;j<=n-1;j++) {
        if (j==i) {continue;}
        temp = temp + a[n*i+j];
    }
    a[n*i+i] = temp+1;
}

//  Solve the system. Break when residue is low
for (k=0;k<=maxiter-1;k++) {
    for (i=0;i<=n-1;i++) {
        temp=0;
        for (j=0;j<=n-1;j++) {
            if (j==i) {continue;}
            temp = temp + a[n*i+j]*x[j];
            }
        xnew[i] = (b[i]-temp)/a[n*i+i];
    }
    temp=0;
    for (i=0;i<=n-1;i++) {
        temp = temp + fabs(x[i]-xnew[i]);
        x[i]=xnew[i];
    }
    if (temp<0.0001) {
        break;
    }
}

printf("Iterations = %d\n",k-1);

return 0;
}

跳出循环的标准非常简单。这个程序永远不应该失败。但它显然没有收敛(它用完了循环中的所有迭代)，除非我将 float 更改为 double 。 float 的精度比这高得多。怎么了？在 Windows 7 下使用 CodeBlocks 16.01 进行编译(如果这很重要的话)。

Answer 1

if (temp<0.0001) {给出的请求太细了float和值。

尝试了不同的限制方法，添加 ULP x[i] 的差异和xnew[i] .

#include <assert.h>
#include <stdint.h>

static uint32_t ULPf(float x) {
  union {
    float f;
    uint32_t u32;
  } u;
  assert(sizeof(float) == sizeof(uint32_t));
  u.f = x;
  if (u.u32 & 0x80000000) {
    u.u32 ^=  0x80000000;
    return    0x80000000 - u.u32;
  }
  return u.u32 + 0x80000000;
}

static uint32_t ULP_diff(float x, float y) {
  uint32_t ullx = ULPf(x);
  uint32_t ully = ULPf(y);
  if (x > y) return ullx - ully;
  return ully - ullx;
}

...

  uint64_t sum0 = -1;
  unsigned increase = 0;
  for (k = 0; k <= maxiter - 1; k++) {
    ...
    uint64_t sum = 0;
    for (i = 0; i <= n - 1; i++) {
      uint32_t e = ULP_diff(x[i], xnew[i]);
      // printf("%u %e %e %llu\n", i, x[i],  xnew[i], (unsigned long long) e);
      sum += e;
      x[i] = xnew[i];
    }
    if (sum < sum0) {
      // answer is converging
      sum0 = sum;
      increase = 0;
    } else {
      increase++;
      // If failed to find a better answer in `n` iterations and 
      //   code did at least n*N iterations, break.
      if (increase > n && k > n*n) break;
    }