分布式稀疏矩阵生成小程序

       生成稀疏矩阵，并统计每行非零元素

#include <stdio.h>
  2 #include <iostream>
  3 
  4 int main()
  5 {
  6         int NumberOfNonzeros = 0;
  7         int totalNumberOfNonzeros = 0;
  8         int nx = 4;
  9         int ny = 4;
 10         int nnzPerRow = 9 ;
 11 
 12         int gnx = nx * ny;
 13  

 14         double** matrixValues = new double*[gnx];
 15         int** mtxIndG = new int*[gnx];
 16         int* nonzerosInRow = new int[gnx];
 17 
 18 /*       initial matrixvalue    */
 19         for(int i=0; i<gnx; i++)
 20         {
 21                 matrixValues[i] = 0;
 22                 mtxIndG[i] = 0;
 23         }
 24 
 25         for(int i=0; i<gnx; i++)
 26         {
 27                 matrixValues[i] = new double[nnzPerRow];
 28                 mtxIndG[i] = new int[nnzPerRow];
 29         }
 30 
 31         /* fill value --> sparse matrix */
 32         for(int iy=0; iy<ny; iy++)
 33         {
 34                 for(int ix=0; ix<nx; ix++)
 35                 {
 36                         int currentRow = iy * nx + ix;
 37 
 38                         char NumberOfNonzerosInRow = 0;
 39                         double* currentValuePointer = matrixValues[currentRow];
 40                         int* currentIndexP = mtxIndG[currentRow];
 <pre name="code" class="cpp">42                         for(int sy=-1; sy<=1; sy++){
 43                                 if(iy+sy>-1 && iy+sy<ny){
 44                                         for(int sx=-1; sx<=1;sx++){
 45                                                 if(ix+sx>-1 && ix+sx<nx){
 46                                                         int curcol = currentRow + sy*nx + sx;
 47                                                         if(curcol == currentRow)
 48                                                         {
 49                                                                 *(currentValuePointer++) = 4;
 50                                                         }else{
 51                                                                 *(currentValuePointer++) = -1.    0;
 52                                                         }
 53                                                         NumberOfNonzerosInRow++;
 54                                                         *currentIndexP++ = curcol;
 55 
 56                                                 }
 57 
 58                                         }
 59                                 }
 60                         }
 61 
 62                         nonzerosInRow[currentRow] = NumberOfNonzerosInRow;
 63                         totalNumberOfNonzeros += NumberOfNonzerosInRow;
 64 
 65                 }
 66         }
 67 
 68 
 69         // output matrix 
 70         FILE* fm = 0;
 71         fm = fopen("matrix.dat", "w");
 72         for(int i=0; i< gnx; i++)
 73         {
 74                 const double* const currentRowValues = matrixValues[i];
 75                 const int* const currentRowIndices = mtxIndG[i];
 76                 const int currentNumberOfNonzeros = nonzerosInRow[i];
 77 
 78                 for(int j=0; j<currentNumberOfNonzeros; j++)
 79                         fprintf(fm, "%d %d %22.16e\n", i+1, currentRowIndices[j]+1, currentRow    Values[j]);
 80         }
 81 
 82         fclose(fm);

 84         // see full matrix in Matlab
 85         // load(matrix.dat)
 86         // spconvert(matrix)
 87 
 88         return 0;
 89 }