Внедрение предварительно обусловленного градиента конъюгата с помощью PETSc

Question

Я пишу код на C для запуска OLS на большом наборе данных, состоящем примерно из 40 миллионов наблюдений. Я заинтересован в использовании предварительно обусловленного алгоритма сопряженного градиента (PCG), использующего неполное разложение Холецкого (ICC) для получения предварительного кондиционера. У меня есть код MATLAB, который делает то же самое правильно и должен использоваться в качестве эталона для проверки правильности сценария Си. Короче говоря, я хочу повторить процедуры MATLAB "ichol" и "pcg", используя C.

В настоящее время я использую библиотеку PETSc C для PCG и ICC. Я использую R для считывания набора данных и вычисления X'X и X'Y, и вывожу оба в текстовые файлы для использования C. Матрица имеет формат разреженных координат (3 вектора - один с индексами строк, один с индексами столбцов и один со значением для всех ненулевых элементов матрицы). Полная матрица имеет размер примерно 8 миллионов на 8 миллионов с около 166 миллионами ненулевых записей.

Я читаю файлы на C, строю матрицу и вектор и запускаю PCG, используя процедуру PETSc KSPSolve, выводя вектор линейных коэффициентов b в текстовый файл. Код завершается без проблем, но решения, которые я получаю, бессмысленны - коэффициенты порядка 2 миллионов, когда я должен получить что-то вроде 2 или 3 (от запуска кода MATLAB).

Однако, когда я ограничиваю набор данных меньшим подвыборкой (например, первые 300000 наблюдений), результаты кода C и кода MATLAB очень хорошо совпадают. После некоторых проб и ошибок я установил, что результаты совпадают, пока я не достигну очень определенного размера выборки (363 855), после чего результаты кода C начинают расходиться. Расхождение не является внезапным взрывом - сначала максимальная абсолютная разница между кодом C и кодом MATLAB составляет порядка 700, но быстро возрастает до 2 миллионов, упомянутых выше. Я посмотрел на внутреннюю статистику мониторинга сходимости PETSc, и сама программа, похоже, считает, что норма истинного остатка (X'Y-X'Xb) снижается до довольно низких уровней (порядка 10^(-3)), а если Я вычисляю остаток вручную в MATLAB, используя вывод из C, получаю, что неудивительно, очень большие числа. Более того, если я возьму подвыборку из последних (а не из первых) 363 855 или даже 1 миллиона наблюдений, то код C и код MATLAB по-прежнему достаточно хорошо совпадают.

Я ищу предложения о том, что может быть не так и какие шаги предпринять, чтобы диагностировать / исправить это. Я довольно плохо знаком с C, поэтому основные ошибки, безусловно, возможны.

Соответствующий раздел кода вставлен ниже:

/* Read in matrix, rhs, and IA, JA, RA vectors from ascii files */

   ierr   = PetscPrintf(PETSC_COMM_SELF,"\n Read matrix in ascii format ...\n");CHKERRQ(ierr);
   ierr   = PetscFOpen(PETSC_COMM_SELF,Ain,"r",&Afile);CHKERRQ(ierr);
   nsizes = 3;
   fscanf(Afile,"%d,%d,%d\n",&m,&n,&nz);
   ierr = PetscPrintf(PETSC_COMM_SELF,"m: %d, n: %d, nz: %d \n", m,n,nz);CHKERRQ(ierr);

/*Allocating matrix and vectors*/

   ierr = MatCreate(PETSC_COMM_SELF,&A);CHKERRQ(ierr);
   ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
   ierr = MatSeqAIJSetPreallocation(A,nz/m,nnz_array);CHKERRQ(ierr);
   ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);


/*Building matrix*/

    for (i=0; i<nz; i++) {
      fscanf(Afile,"%g,%g,%le\n",&_row,&_col,(double*)&val);
      row = (PetscInt)_row;
      col = (PetscInt)_col;
      ierr = MatSetValues(A,1,&row,1,&col,&val,INSERT_VALUES);
      CHKERRQ(ierr);

    }

    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    fflush(stdout);
    fclose(Afile);
    free(nnz_array);

  /*Reading-in X'Y vector*/
  ierr = VecCreate(PETSC_COMM_SELF,&b);CHKERRQ(ierr);
  ierr = VecSetSizes(b,PETSC_DECIDE,n);CHKERRQ(ierr);
  ierr = VecSetFromOptions(b);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_SELF,"\n Read rhs in ascii format ...\n");CHKERRQ(ierr);
  ierr = PetscFOpen(PETSC_COMM_SELF,rhs,"r",&bfile);CHKERRQ(ierr);
    for (i=0; i<n; i++) {
      fscanf(bfile,"%d,%le\n",&dummy,(double*)&val);
      ierr = VecSetValues(b,1,&i,&val,INSERT_VALUES);CHKERRQ(ierr);
    }
    ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
    fflush(stdout);
    fclose(bfile);

  /*Creating vector to store solution in*/
  ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) x, "Solution");CHKERRQ(ierr);
  ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                Create the linear solver and set various options
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /*
     Create linear solver context
  */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);

  /*
     Set operators. Here the matrix that defines the linear system
     also serves as the preconditioning matrix.
  */
  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);

  /*
     Set linear solver defaults for this problem.
  */
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCICC);CHKERRQ(ierr);
  ierr = KSPSetTolerances(ksp,1.e-10,1.e-10,PETSC_DEFAULT,1000);CHKERRQ(ierr);
  ierr = KSPSetType(ksp,KSPCG);
  /*
    Set runtime options.
  */
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                      Solve the linear system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /*
     Solve linear system
  */
    printf("\nBeginning to solve system:\n"); 
    ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);


  /*
     View solver info.
  */
  ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);