0

我正在编写一些 C 代码来在大约 4000 万个观测值的大型数据集上运行 OLS。我对使用不完全 Cholesky (ICC) 分解来获得预处理器的预处理共轭梯度算法 (PCG) 很感兴趣。我有一个 MATLAB 代码可以正确执行相同的操作，并将其用作基准以确保 C 脚本正确。简而言之，我想使用 C 复制 MATLAB 的“ichol”和“pcg”程序。

目前我正在为 PCG 和 ICC 使用 PETSc C 库。我使用 R 读入数据集并计算 X'X 和 X'Y，并将两者都输出到文本文件以供 C 使用。矩阵采用坐标稀疏格式（3 个向量 - 一个带有行索引，一个带有列索引，和一个具有所有非零矩阵条目的值）。完整矩阵的大小约为 800 万乘 800 万，其中包含约 1.66 亿个非零条目。

我在 C 中读入文件，构造矩阵和向量，并使用 PETSc 的程序 KSPSolve 运行 PCG，将线性系数 b 的向量输出到文本文件。代码顺利完成，但我得到的解决方案是荒谬的 - 系数约为 200 万，而我应该得到 2 或 3 之类的东西（通过运行 MATLAB 代码）。

但是，当我将数据集限制为较小的子样本（例如，前 300,000 个观察值）时，C 代码和 MATLAB 代码的结果非常吻合。经过反复试验，我确定结果一致，直到达到非常特定的样本大小 (363,855)，此时 C 代码结果开始出现差异。分歧并不是突然爆发的——起初 C 代码和 MATLAB 代码之间的最大绝对差值在 700 左右，但很快增长到上面提到的 200 万。我查看了 PETSc 的内部收敛监控统计数据，程序本身似乎认为真实残差 (X'YX'Xb) 的范数下降到相当低的水平（10^(-3) 的顺序），而如果我计算使用 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);
4

0 回答 0