Malloc ошибка освобождения памяти - C (СИ)

double SizeArr = 5, wData = 8, **Arr0, *Arr1, *Arr2;
 
Len = SizeArr * sizeof(double);
Arr0 = (double**)malloc(SizeArr * sizeof(*Arr0));
    for (i = 0; i < SizeArr; ++i) Arr[i] = (double*)malloc(Len);
 
Arr1 = (double*)malloc(Len);
Arr2 = (double*)malloc(Len);
 
for (i = 0; i < SizeArr; ++i) 
    Arr[i] = (double*)realloc(Arr[i], wData * sizeof(*Arr));
 
free(Arr1); Arr1 = NULL; // нет ошбки
free(Arr2); Arr2 = NULL; // ошибка
for (i = 0; i < SizeArr; ++i) 
    free(Arr0[i]); // ошибка - срабатывание точки останова
free(Arr0); Arr0 = NULL;

#include <math.h>
#include <stdio.h>
#include <iostream>
#include <ctime>
 
void mxmt(double **Data, double **Arr, size_t wData, size_t SizeArr);
void tred2(double **Arr, size_t SizeArr, double *D0, double *D1);
void tqli1(double **Arr, size_t SizeArr, double *D0, double *D1, size_t *MaxIter);
void bevm(double **Arr, size_t SizeArr, double *D0, double *D1, double **DT);
void ivm(double **DT, double *D0, size_t SizeArr);
void vxb(double *A, double **B, double **C, size_t whM, size_t wB);
void smdd(double **C, double *D0, size_t SizeArr);
void axb(double **A, double **B, double **C, size_t hA, size_t whM, size_t wB);
void output1(double *D, size_t SizeArr);
void output2i(double **Arr, size_t SizeArr);
void output2t(double **Arr, size_t SizeArr);
void output2r(double **Arr, size_t SizeArr, size_t wData);
void output2rt(double **Arr, size_t SizeArr, size_t wData);
 
int main() {
    setlocale(0, "ru-ru"); srand((size_t)time(0));
    double **Data, **Arr, **DT, *D0, *D1, **SU;
    size_t i, SizeArr = 5, wData = 8, Len, *MaxIter = &i;
  
    Len = wData * sizeof(**Data);
    Data = (double**)malloc(SizeArr * sizeof(*Data));
        for (i = 0; i < SizeArr; ++i) Data[i] = (double*)malloc(Len);
    
    for (i = 0; i < SizeArr; ++i) 
        for (size_t j = 0; j < wData; ++j) Data[i][j] = 5 - rand() % 11;
 
    Len = SizeArr * sizeof(double);
    Arr = (double**)malloc(SizeArr * sizeof(*Arr));
        for (i = 0; i < SizeArr; ++i) Arr[i] = (double*)malloc(Len);
 
    mxmt(Data, Arr, wData, SizeArr); 
    puts("Исходные данные (Data): "); output2r(Data, SizeArr, wData);
    puts("Исходная левая симметричная матрица (Data x Data^T): "); output2i(Arr, SizeArr);
 
    D0 = (double*)malloc(Len); //memset(&D0[0], 0, Len);
  D1 = (double*)malloc(Len); //memset(&D1[0], 0, Len);
    DT = (double**)malloc(3 * sizeof(**DT));
        for (i = 0; i < 3; i++) DT[i] = (double*)malloc(Len);
 
    tred2(Arr, SizeArr, D0, D1);
    for (i = 0; i < SizeArr; ++i) { DT[0][i] = D0[i]; DT[1][i] = D1[i]; }
    i = 30; tqli1(Arr, SizeArr, D0, D1, MaxIter);
    for (i = 0; i < SizeArr; ++i) DT[2][i] = D0[i];
    puts("Собственные значения (L): "); output1(D0, SizeArr);
    bevm(Arr, SizeArr, D0, D1, DT);
    puts("Собственные векторы (V): "); output2t(Arr, SizeArr);
    ivm(DT, D0, SizeArr);
    puts("Обратный вектор сингулярных значений (S'): "); output1(D0, SizeArr);
 
    // Arr = D0 x Arr x Data [ V^(T) = (S^(T))^(-1) x U^(T) x A ]
    SU = (double**)malloc(SizeArr * sizeof(**SU));
        for (i = 0; i < SizeArr; i++) SU[i] = (double*)malloc(Len);
    
    // используем формулу умножения вектора на не транспонированную матрицу, 
    // поскольку её содежимое уже транспонировано
    vxb(D0, Arr, SU, SizeArr, SizeArr); //output2t(SU, SizeArr);
    smdd(SU, D0, SizeArr); //output2t(SU, SizeArr);
 
    for (i = 0; i < SizeArr; ++i) 
        Arr[i] = (double*)realloc(Arr[i], wData * sizeof(*Arr));
 
    axb(SU, Data, Arr, SizeArr, SizeArr, wData);
    puts("Матрица векторов (V), нулевые вектора опущены: "); output2rt(Arr, wData, SizeArr);
 
    system("pause");
 
    for (i = 0; i < SizeArr; ++i) free(Data[i]); free(Data); Data = NULL;
    for (i = 0; i < SizeArr; ++i)   free(SU[i]); free(SU); SU = NULL;
    for (i = 0; i < 3; ++i) free(DT[i]); free(DT); DT = NULL;
    free(D0); D0 = NULL;
    //free(D1); //D1 = NULL;
    for (i = 0; i < SizeArr; ++i) 
        free(Arr[i]); 
    //free(Arr); Arr = NULL;
    
    return 0;
}
 
void mxmt(double **Data, double **Arr, size_t wData, size_t SizeArr) {
    for (size_t i = 0; i < SizeArr; ++i)
        for (size_t j = i; j < SizeArr; ++j) {
            Arr[j][i] = 0.;
            for (size_t k = 0; k < wData; ++k)
                Arr[j][i] += Data[i][k] * Data[j][k];
            Arr[i][j] = Arr[j][i];
        }
}
 
void tred2(double **Arr, size_t SizeArr, double *D0, double *D1) {
    size_t i, j, k, l; double Scale, Sum, Pre, HH;
    
    /* получение ортогональной матрицы трансформации */
    /* проход по стадиям/итерациям процесса редукции */
    for (i = SizeArr - 1; i > 0; --i) {
      l = i - 1; D0[i] = Scale = 0.;
    
      /* вычислить сумму модулей строки под главной диагональю */
      for (j = 0; j < i; ++j) Scale += fabs(Arr[i][j]);
    
      /* сложный процесс везде, кроме последней стадии/итерациям, весь 
      обрабатываемый диапазон находится под главной диагональю матрицы */
      if (i < 2 || !Scale) D1[i] = Arr[i][l];
      else { /* преобразование матрицы */
        /* масштабируем строку */
        for (j = 0; j < i; ++j) Arr[i][j] /= Scale;
    
        /* вычислить сумму квадратов строки */
        for (j = 0; j < i; ++j) D0[i] += Arr[i][j] * Arr[i][j];
    
        /* работаем с вектором U */
        D1[i] = (Arr[i][l] < 0. ? sqrt(D0[i]) : -sqrt(D0[i]));
        D0[i] -= Arr[i][l] * D1[i]; Arr[i][l] -= D1[i];
        D1[i] *= Scale;
    
        /* вычисляем предварительный коэффициент Pre */
        for (Pre = 0., j = 0; j < i; ++j) {
          Arr[j][i] = Arr[i][j] / D0[i];
          for (Sum = 0., k = 0; k <= j; ++k) Sum += Arr[j][k] * Arr[i][k];
          for (k = j + 1; k < i; ++k) Sum += Arr[k][j] * Arr[i][k];
          D1[j] = Sum / D0[i]; Pre += Arr[i][j] * D1[j];
        }
        /* трансформируем матрицу */
        HH = Pre / (D0[i] + D0[i]);
        for (j = 0; j < i; ++j) {
          D1[j] -= HH * Arr[i][j];
          for (k = 0; k <= j; ++k)
            Arr[j][k] -= (D1[k] * Arr[i][j] + D1[j] * Arr[i][k]);
        }
      }
    }
    /* сохраняем главную диагональ */
    for (D1[0] = 0., i = 0; i < SizeArr; ++i) D0[i] = Arr[i][i]; /**/
}
 
void tqli1(double **Arr, size_t SizeArr, double *D0, double *D1, size_t *MaxIter) {
  size_t i, j, k, Iter; double b, c, DD, f, g, p, r, s;
 
  /* сдвигаем массив влево для удобства */
  for (i = 1; i < SizeArr; ++i) D1[i - 1] = D1[i]; D1[i - 1] = 0.;
      
  for (i = 0; i < SizeArr; ++i) {
    Iter = *MaxIter;
    do { /* диагонализируем миноры с текущей строки */
      /* находим малый поддиагональный элемент для расщепления */
      for (j = i; j < SizeArr - 1; ++j) {
        DD = fabs(D0[j]) + fabs(D0[j + 1]);
        if ((double)(fabs(D1[j] + DD)) == DD) break;
      }
      /* если верхний левый угол минора еще не диагональный */
      if (j != i) {
        /* формирование сдвига */
        g = (D0[i + 1] - D0[i]) / (2. * D1[i]); r = hypot(1., g);
        /* находим сумму модулей d_m и k_s с сохранением исходного знака */
        g += (g >= 0. ? fabs(r) : -fabs(r)); g = D0[j] - D0[i] + D1[i] / g;
        /* плоская ротация оригинального QL алгоритма, сопровождаемая
        ротациями Гивенса для восстановления трехдиагональной формы */
                for (s = c = 1., p = 0., k = j - 1; k < j && k >= i; --k) {
          f = s * D1[k]; b = c * D1[k]; D1[k + 1] = r = hypot(f, g);
          if (r == 0.) { D0[k + 1] -= p; D1[j] = 0.; break; }
          /* основные действия на ротации */
          s = f / r; c = g / r; g = D0[k + 1] - p;
          r = (D0[k] - g) * s + 2. * c * b;
          p = s * r; D0[k + 1] = g + p; g = c * r - b;
        }
        /* безусловный переход к новой итерации при нулевом знаменателе
        и недоведенной до конца последовательности ротаций */
        if (r == 0. && k < j && k >= i) continue;
        /* сохраняем новые значения диагоналей */
        D0[i] -= p; D1[i] = g; D1[j] = 0.;
      }
    } while (j != i && --Iter);
    if (!Iter) { printf("Привышен заданный порог итераций на цикле %d \n", i); break; }
  }
}
 
void bevm(double **Arr, size_t SizeArr, double *D0, double *D1, double **DT) {
    double Div, Delta; size_t i, j, k, l, Rang;
    // цикл поиска сбственных векторов для каждого собственого значения
    for (i = 0; i < SizeArr; ++i) {
        // получаем главную и прилигающие к ней диагонали
        for (j = 0; j < SizeArr; ++j) { D0[j] = DT[0][j] - DT[2][i]; D1[j - 1] = DT[1][j]; }
        // получаем триугольную/трапциидальную матрицу из трёхдиагональной
        Div = 1.; Rang = l = SizeArr - 1; D1[l] = D1[0];
        for (j = 0, k = 1; k < l; ++j, ++k) {
            D0[k] = (D0[j] * D0[k] - D1[l] * D1[j]) / Div;
            D1[l] = D1[k]; D1[k] = D0[j] / Div * D1[k]; Div = D0[j];
        }
        // подсчитываем среднюю линейную погрешность с некоторым запасом
        // кого это не устраивает - считает погрешность как положено
        D0[k] = (D0[j] * D0[k] - D1[l] * D1[j]) / Div;
        Delta = 1.00059 * D0[k]; Div = D0[k] / k; D0[k] = 1.;
        // поиск свободных значений / подсчёт ранга матрицы
        for (j = SizeArr - 2; j; --j, Delta -= Div)
            if (fabs(D0[j]) < Delta && fabs(D1[j]) < Delta) { D0[j] = 1.; --Rang; }
            else break;
        // получаем собственный вектор для текущего собственного значения
        for (j = Rang - 1, k = Rang; j < Rang; --j, --k) D0[j] = -D0[k] / D0[j] * D1[j];
        // нормализуем вектор
        Div = 0.; for (j = 0; j < SizeArr; ++j) Div += D0[j] * D0[j];
        Div = sqrt(Div); for (j = 0; j < SizeArr; ++j) D0[j] /= Div;
        // сохраняем результат в массив векторов в транспонированном виде (по строкам)
        for (j = 0; j < SizeArr; ++j) Arr[i][j] = D0[j];
    }
}
 
void ivm(double **DT, double *D0, size_t SizeArr) {
    size_t i; // подсчитываем сингулярные значения
    for (i = 0; i < SizeArr; ++i) D0[i] = (DT[2][i] > 0 ? sqrt(DT[2][i]) : 0);
    // находим обратную вектор матрицу
    for (i = 0; i < SizeArr; ++i) D0[i] = 1. / D0[i];
}
 
void vxb(double *V, double **B, double **C, size_t whM, size_t wB) {
    size_t i, j;
    for (i = 0; i < whM; ++i) for (j = 0; j < wB; ++j) C[i][j] = V[i] * B[i][j];
}
 
void smdd(double **C, double *D0, size_t SizeArr) {
    // сортировка строк матрицы по убыванию значений вектор матрицы
    size_t Pos = 0, i, j; double Max, Swp, *pSwp;
    for (i = 0; i < SizeArr - 1; ++i) {
        for (j = i + 1, Max = D0[i]; j < SizeArr; ++j)
            if (Max < D0[j]) Max = D0[j], Pos = j;
        if (i != Pos) {
            Swp = D0[i]; D0[i] = D0[Pos]; D0[Pos] = Swp;
            pSwp = C[i]; C[i] = C[Pos]; C[Pos] = pSwp;
        }
    }
}
 
void axb(double **A, double **B, double **C, size_t hA, size_t whM, size_t wB) {
    for (size_t i = 0; i < hA; ++i)
        for (size_t j = 0; j < wB; ++j) {
            C[i][j] = 0.; 
            for (size_t k = 0; k < whM; ++k) C[i][j] += A[i][k] * B[k][j];
        }
}
 
void output1(double *a, size_t n) {
  for (size_t x = 0; x < n; ++x) printf("%14.10lf ", a[x]); puts("\n");
}
 
void output2i(double **a, size_t n) {
    for (size_t y = 0; y < n; ++y) {
        for (size_t x = 0; x < n; ++x) printf("%5.0lf ", a[y][x]); puts("");
    } puts("");
}
 
void output2t(double **a, size_t n) {
  for (size_t y = 0; y < n; ++y) {
    for (size_t x = 0; x < n; ++x) printf("%14.10lf ", a[x][y]); puts("");
  } puts("");
}
 
void output2r(double **a, size_t n, size_t m) {
    for (size_t y = 0; y < n; ++y) {
        for (size_t x = 0; x < m; ++x) printf("%5.0lf ", a[y][x]); puts("");
    } puts("");
}
 
void output2rt(double **a, size_t n, size_t m) {
    for (size_t y = 0; y < n; ++y) {
        for (size_t x = 0; x < m; ++x) printf("%14.10lf ", a[x][y]); puts("");
    } puts("");
}