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

#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 }