Jacobi 方法适用于 double，不适用于 float。怎么了？

Question

我编写了一个小程序来使用 Jacobi（迭代）方法求解一个包含 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;
}

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

Answer 1

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

通过添加 x[i] 和 xnew[i] 之差的 ULP 尝试了不同的限制方法。

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

Answer 2

似乎 float 数据类型没有上述算法所需的精度，鉴于编码方式。该算法确实会收敛，但是 "residue" 永远不会低到足以退出循环。

我的理解是，由于 float 变量在内部存储的方式，你不能用极小的 (0.0001) 和极大的 (RAND_MAX) 进行计算数字并期望合理的准确性，如上例所示（temp 在最内层循环中增长到一个巨大的数字）。

因此，设置b[i] = rand()%60;和x[i] = rand()%60;将缓解问题。

设置 b[i] = rand()%6; x[i] = rand()%6; 和 a[n*i+j] = rand()%6 将最终满足更严格的退出循环条件。

有趣的是，建立更大的对角线优势（将 a[n*i+i] = temp+1' 更改为 a[n*i+i] = temp+10; 也会使程序收敛，而以前不会收敛。

我不熟悉其他人描述的 ULP 条件，但会花一些时间了解它

如果未来的读者有时间和精力，也许他们应该阅读 "What Every Computer Scientist Should Know About Floating-Point Arithmetic"，即使我没有。

顺便说一句，xorg 是存储原始 x 向量，用于调试目的，因为我很难让 CodeBlocks 进行调试

感谢大家的贡献。

Jacobi 方法适用于 double，不适用于 float。怎么了？

Jacobi method works with double, fails with float. What is wrong?

c

precision

double