SVD MPI 实现

mzg12345678

于 2018-01-03 16:23:42 发布

阅读量295

点赞数

run_lsi.sh
1 ROOT_PATH=/da1/search/hpc_install/install
2 export LD_LIBRARY_PATH=$ROOT_PATH/mpich/lib/:$LD_LIBRARY_PATH
3 export PATH=$ROOT_PATH/mpich/bin:$PATH
4 export LD_LIBRARY_PATH=/da1/search/gcc-4.8/lib64:/da1/search/gcc-4.8/lib:$LD_LIBRARY_PATH
5
6
7 mkdir -p error_out
8 mkdir -p stdout
9
10 pscp -h host /da1/search/lsi_workspace/libcrypto.so.6 /da1/search/lsi/doc2term/bin
11 pscp -h host /da1/search/lsi_workspace/libquadmath.so.0 /da1/search/lsi/doc2term/bin
12 pscp -h host -t 0 -e error_out -o stdout $1 "/da1/search/lsi/doc2term/bin"
13 mpirun -f host -n 90 -wdir /da1/search/lsi/doc2term/bin ./$1 -svd_nsv 100 -svd_ncv 200 -svd_max_it 200 -svd_lanczos_oneside

Petsc + slepc

2G 的网页的title 和 qanchor 数据，2M 的term，大概有30G个非0的数据

// 每台机器读取自己part的数据，加载进内存，并计算整个矩阵的大小

  RowBasedMatrix rbm;
  init_row_based_matrix(&rbm);
  PetscInt local_max_col = 0;
  PetscInt local_row = load_matrix(&rbm, &local_max_col);
  assert(rbm.total_row_num == local_row);

  PetscScalar* array_value = (PetscScalar*)malloc(rbm.max_col_num * sizeof(PetscScalar));
  for (i = 0; i < rbm.max_col_num; i++) {
    array_value[i] = 1.0;
  }

  ierr = PetscPrintf(PETSC_COMM_WORLD,"load matrix from file finished, max_col\n");CHKERRQ(ierr);
  int max_col_index = 0;
  MPI_Status status;

  PetscInt local_row_start = 0;

  for (i = rank + 1; i < num_procs; i ++) {
    MPI_Send(&local_row, 1, MPI_UNSIGNED_LONG, i, 3, PETSC_COMM_WORLD);
  }

  if (rank != 0) {
    MPI_Send(&local_max_col, 1, MPI_INT, 0, 1, PETSC_COMM_WORLD);
    PetscInt other_row_start = 0;
    for (i = 0; i < rank; i++) {
        MPI_Recv(&other_row_start, 1, MPI_INT, i, 3, PETSC_COMM_WORLD, &status);
        local_row_start += other_row_start;
    }
    MPI_Recv(&max_col_index, 1, MPI_INT, 0, 2, PETSC_COMM_WORLD, &status);
  } else {
      int other_col_num = 0;
      max_col_index = local_max_col;
      for (i = 1; i < num_procs; i++) {
        MPI_Recv(&other_col_num, 1, MPI_INT, i, 1, PETSC_COMM_WORLD, &status);
        if (max_col_index < other_col_num) {
            max_col_index = other_col_num;
        }
      }

      for (i = 1; i < num_procs; i++) {
        MPI_Send(&max_col_index, 1, MPI_INT, i, 2, PETSC_COMM_WORLD);
      }
  }
 // 变为矩阵的结构
  calc_dnz_onz(&rbm, local_row, global_col_start, global_col_end, dnz, onz);
  ierr = MatCreateAIJ(PETSC_COMM_WORLD, local_row, PETSC_DECIDE, PETSC_DETERMINE, max_col_index + 1, 0, dnz, 0, onz, &A);
  //MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
  //ierr = MatCreate(PETSC_COMM_WORLD, &A);CHKERRQ(ierr);
  //ierr = MatSetSizes(A, local_row, local_row, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);

  //ierr = MatSetUp(A); CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD, "total_element_num:%d, total_row_num:%d\n", rbm.total_element_num, rbm.total_row_num);
  //ierr = MatMPIAIJSetPreallocation(A, 3 * rbm.total_element_num / rbm.total_row_num, PETSC_NULL, 3 * rbm.total_element_num / rbm.total_row_num, PETSC_NULL);
  //ierr = MatMPIAIJSetPreallocation(A, 0, dnz, 0, onz);

  MatGetOwnershipRange(A, &global_row_start, &global_row_end);

  ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "rank=%d, row:%d, col:%d, start:%d, global_start:%d, global_row_num:%d, gcs:%d, gce:%d\n",
            rank, local_row, local_max_col, local_row_start, global_row_start, global_row_end - global_row_start, global_col_start, global_col_end);

  PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout);

  //MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);

  CHKERRQ(ierr);

  PetscPrintf(PETSC_COMM_WORLD, "row %d, col %d,first:%llu, last:%llu\n",
        0, rbm.total_row_data[0].col_num, rbm.total_row_data[0].col_ids[0],
        rbm.total_row_data[0].col_ids[rbm.total_row_data[0].col_num - 1]);

  for (i = 0; i < global_row_end - global_row_start; i++) {
    if (rbm.total_row_data[i].col_num == 0) {
        continue;
    }
    //MatSetValues(A, 1, &row_id, head_data->row_num, head_data->row_ids, col_data, INSERT_VALUES);
    PetscInt target_global_row = i + global_row_start;
    //MatSetValues(A, 1, &global_row, rbm.total_row_data[i].col_num, rbm.total_row_data[i].col_ids, col_data, INSERT_VALUES);
    MatSetValues(A, 1, &target_global_row, rbm.total_row_data[i].col_num, rbm.total_row_data[i].col_ids, array_value, INSERT_VALUES);
    free(rbm.total_row_data[i].col_ids);
    //free(rbm.total_row_data[i].col_values);

    assert(target_global_row >= global_row_start);
    assert(target_global_row < global_row_end);
  }

 // 计算svd
   ierr = SVDCreate(PETSC_COMM_WORLD,&svd);CHKERRQ(ierr);

  /*
     Set operator
  */
  ierr = SVDSetOperator(svd, A);CHKERRQ(ierr);

  /*
     Set solver parameters at runtime
  */
  ierr = SVDSetFromOptions(svd);CHKERRQ(ierr);
  SVDSetType(svd, SVDLANCZOS);
  SVDLanczosSetOneSide(svd, PETSC_TRUE);
  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                      Solve the singular value system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  BV U, V;
  ierr = BVCreate(PETSC_COMM_WORLD, &U); CHKERRQ(ierr);
  ierr = BVCreate(PETSC_COMM_WORLD, &V); CHKERRQ(ierr);
  ierr = SVDSetBV(svd, U, V); CHKERRQ(ierr);
  ierr = SVDSetImplicitTranspose(svd, PETSC_TRUE); CHKERRQ(ierr);
  ierr = SVDSolve(svd);CHKERRQ(ierr);

  gettimeofday(&t1, NULL);
  PetscInt svd_time = time_diff(&t2, &t1);
  ierr = PetscPrintf(PETSC_COMM_WORLD," svd sovle finished %llu us\n",svd_time);CHKERRQ(ierr);

  ierr = SVDGetIterationNumber(svd,&its);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);CHKERRQ(ierr);

  /*
     Optional: Get some information from the solver and display it
  */
  ierr = SVDGetType(svd,&type);CHKERRQ(ierr);