C++ createMatrix函数代码示例

// calculates \sum_K v_i'*A*v_i
double dosum(Matrix A, Matrix v)
{
  double alpha=0;
  int t = omp_get_max_threads();
  Matrix temp = createMatrix(A->rows,t);
#pragma omp parallel for schedule(static) reduction(+:alpha)
  for( int i=0;i<v->cols;++i ) {
    myMxV(temp->col[omp_get_thread_num()],A,v->col[i]);
    alpha += myinnerproduct(temp->col[omp_get_thread_num()],v->col[i]);
  }
  freeMatrix(temp);

  return alpha;
}

int readTwoOnDimMatrixFromAFile(char *fileName , oneDimMatrix * a , oneDimMatrix * b){
   FILE *fp;
   int numberOfLines = 0;

   fp = openAFileForReading(fileName);
   numberOfLines = readMatrixFromFile(fp, numberOfLines, DETECT_SIZE, a);
   printf("%d  %d %d \n" ,numberOfLines ,  a->nrows , a->ncols );
   numberOfLines= readMatrixFromFile(fp, numberOfLines, DETECT_SIZE , b);
   printf("%d  %d %d \n" ,numberOfLines , b->nrows , b->ncols);
   if(a->nrows != b->nrows)
      inconsistentErrorReport("Different number of rows in each matrix.", 
                              numberOfLines, a->nrows , a->ncols);
   fclose(fp);

   *a = createMatrix(a->nrows , a->ncols);
   *b = createMatrix(b->nrows , b->ncols);
   fp = openAFileForReading(fileName);
   numberOfLines = readMatrixFromFile(fp, numberOfLines, READ_MATRIX, a);
   numberOfLines= readMatrixFromFile(fp, numberOfLines, READ_MATRIX , b);
   fclose(fp);   
   /*   printOneDimMatrix(a);  printf("\n");  printOneDimMatrix(b); */
   return numberOfLines;
}

// creates an identity matrix with numRow rows and numCol columns; works similar to MATLABs eye()
void eye(struct Matrix *eyeMatrix,int numRow,int numCol) {
  createMatrix(eyeMatrix,numRow,numCol);                // create matrix eyeMatrix
  int i,j;
  for (i=0;i<numRow;i++) {                        // iterate through rows
    for (j=0;j<numCol;j++) {                      // iterate through columns
      if (i==j) {
        eyeMatrix->value[i][j]=1.0;                   // if on the diagonal set to 1
      }
      else {
        eyeMatrix->value[i][j]=0.0;                   // else set value to zero
      }
    }
  }
}

int main(int argc, char** argv)
{
  int rank, size;
  init_app(argc, argv, &rank, &size);

  if (argc < 2) {
    printf("usage: %s <N> [L]\n",argv[0]);
    close_app();
    return 1;
  }

  /* the total number of grid points in each spatial direction is (N+1) */
  /* the total number of degrees-of-freedom in each spatial direction is (N-1) */
  int N  = atoi(argv[1]);
  int M  = N-1;
  double L=1.0;
  if (argc > 2)
    L = atof(argv[2]);

  double h = L/N;
  poisson_info_t ctx;
  ctx.A = createPoisson1D(M);

  Vector grid = createVector(M);
  for (int i=0;i<M;++i)
    grid->data[i] = (i+1)*h;

  Matrix u = createMatrix(M, M);
  evalMesh(u->as_vec, grid, grid, poisson_source);
  scaleVector(u->as_vec, h*h);

  double time = WallTime();
  cg(evaluate, u, 1.e-6, &ctx);

  evalMesh2(u->as_vec, grid, grid, exact_solution, -1.0);
  double max = maxNorm(u->as_vec);

  if (rank == 0) {
    printf("elapsed: %f\n", WallTime()-time);
    printf("max: %f\n", max);
  }

  freeMatrix(u);
  freeVector(grid);
  freeMatrix(ctx.A);

  close_app();
  return 0;
}

void printEigen( Matrix *eigen_values, Matrix *eigen_vectors, int num_of_file ) {
	char fn[256];

	Matrix *eigen = createMatrix( eigen_vectors->height + 1, eigen_values->height );
	for ( int row = 0; row < eigen->height; row++ ) {
		eigen->a[row][0] = eigen_values->a[row][row];
		for ( int col = 1; col < eigen->width; col++ ) {
			eigen->a[row][col] = eigen_vectors->a[col-1][row];
		}
	}

	sort( eigen );
	sprintf( fn, "eigen%02d.txt", num_of_file );
	printMatrixToFile( fn, eigen );
}

static ZSolveMatrix Matrix2zsolve(Matrix *M)
{
    int i, j;
    ZSolveMatrix zmatrix;

    zmatrix = createMatrix(M->NbColumns-2, M->NbRows);
    for (i = 0; i < M->NbRows; ++i)
	for (j = 0; j < M->NbColumns-2; ++j) {
	    assert(mpz_cmp_si(M->p[i][1+j], -MAXINT) > 0);
	    assert(mpz_cmp_si(M->p[i][1+j], MAXINT) < 0);
	    zmatrix->Data[i*zmatrix->Width+j] = mpz_get_si(M->p[i][1+j]);
	}

    return zmatrix;
}

Matrix createPoisson1D(int M)
{
  Matrix result = createMatrix(M, M);
  result->data[0][0] = 2.0;
  result->data[1][0] = -1.0;
  for (int i=1;i<M-1;++i) {
    result->data[i][i] = 2.0;
    result->data[i+1][i] = -1.0;
    result->data[i-1][i] = -1.0;
  }
  result->data[M-1][M-1] = 2.0;
  result->data[M-2][M-1] = -1.0;

  return result;
}

SmithWatermanDP::SmithWatermanDP(const char* s1, size_t n1, const char* s2, size_t n2) 
  // recall that the SW matrix has n+1 and m+1 rows and columns respectively
  : DynamicProgramming(n1+1,n2+1)
{
  this->x = s1;
  this->y = s2;
  // default gap penalty is 1 and efault sim matrix is 1 on
  // diagonal and 0 elsewere
  this->gapPenalty = 1;
  this->sim = createDefaultSimilarityMatrix(256);
  // defualt value for backtracking is not enbaled
  this->btEnabled = false;
  this->btMatrix = NULL;
  createMatrix();
}

int main(int argc, char** argv)
{
  int N, K, i, j;
  Matrix A,v;
  double time, sum;

  if (argc < 3) {
    printf("need two parameters, the matrix size and the number of vectors\n");
    return 1;
  }

  N=atoi(argv[1]);
  K=atoi(argv[2]);

  A = createMatrix(N,N);
  // identity matrix
  for (i=0;i<N;++i)
    A->data[i][i] = 1.0;

  v = createMatrix(N,K);

  // fill with column number
  for (i=0;i<K;++i)
    for (j=0;j<N;++j)
      v->data[i][j] = i;

  time = WallTime();
  sum = dosum(A,v);

  printf("sum: %f\n", sum);
  printf("elapsed: %f\n", WallTime()-time);
  freeMatrix(v);
  freeMatrix(A);

  return 0;
}

int main()
{
    FILE * f = fopen("matrix.txt","r");
    readFromAdjMatrix(f);
    printAdjMatrix();
    NodeT** v=(NodeT**)malloc(nrOfVerteces*sizeof(NodeT));
    v=CreateList(adjMatrix);
    int **matrix=(int**)malloc(sizeof(int*)*nrOfVerteces);
    matrix=createMatrix(v);
    printAllList(v);
    printMatrix(matrix);
    bfs(0);
    dfs(0);
    dfsRecurs(0);
    return 0;
}

// sharpening image - apply laplacian filter and
//subtract the resultant image from the original image
struct matrix* imsharpen(struct matrix* image)
{
	struct matrix* filter;
	filter = createMatrix(3,3);

	initMatrix(filter, (void*)laplacianFilter);

	struct matrix* filteredImage = imfilter(image,filter,CORRELATION_OPERATION,MODE_PAD_ZERO);

	struct matrix* sharpenedImage = subtractMatrices(image,filteredImage);

	destroyMatrix(filter);
	destroyMatrix(filteredImage);

	return sharpenedImage;
}

double HMM::backwardAlgorithm(vector<int> words)
{
	int len = words.size();
	double ** backward = createMatrix(len, _maxState);
	computeBackwardMatrix(words, backward, len);

	double p = 0;
	for (int i = 0; i < _maxState; i++)
	{
		p += backward[0][i] * _pObservation[i][words[0]] * initial_probability[i];
	}

	freeMatrix(backward, len, _maxState);

	return p;
}

CvMatND* FeatPyramid::padArray (CvMatND *mat, int dimPad[3], float val)
{
  int dims[3];

  if (mat == NULL)
  {
    for (int i = 0; i < 2; i++)
      dims[i] = dimPad[i] * 2;

    dims[2] = 32 + (dimPad[2] * 2);
  }

  else
  {
    for (int i = 0; i < 3; i++)
      dims[i] = mat->dim[i].size + (dimPad[i] * 2);
  }

  // New bigger matrix is created
  assert (dims[0] > 0);
  assert (dims[1] > 0);
  assert (dims[2] > 0);

  //CvMatND *padMat = createNDMatrix (3, dims, CV_64FC1);
	CvMatND *padMat;
	createMatrix (3, dims, CV_64FC1, &padMat);
  assert (padMat != NULL);

  // Firstable, all the elements are setted to val
  for ( int i = 0; i < padMat->dim[0].size; i++ )
    for ( int j = 0; j < padMat->dim[1].size; j++ )
      for ( int k = 0; k < padMat->dim[2].size; k++ )
        cvSetReal3D (padMat, i, j, k, val);

  if (mat != NULL)
  {
    // Original values are setted
    for ( int i = dimPad[0]; i < padMat->dim[0].size - dimPad[0]; i++ )
      for ( int j = dimPad[1]; j < padMat->dim[1].size - dimPad[1]; j++ )
        for ( int k = dimPad[2]; k < padMat->dim[2].size - dimPad[2]; k++ )
          cvSetReal3D (padMat, i, j, k, cvGetReal3D (mat, i-dimPad[0], j-dimPad[1], k-dimPad[2]));
  }

  return padMat;
}

NS_CORE InternalGearRatios DefKNuton::solveNuton( const Jacobi& jacobian, System& system )
{
	const double eps = 0.001;
	const int maxIterCount = 100;
	double norm;
	int iterCount = 0;
	bool notFinded = false;

	do
	{
		// создаем матрицу с уравнениями
		auto matrix = createMatrix( jacobian, system );
		// создаем вектор правых частей уравнений
		auto rightParts = createRightParts( system );
		// решаем систему уравнений и получаем дельту
		auto next = MatrixOperations::solveGaus( matrix, rightParts );
		if ( next.size() == 0 )
		{
			next = MatrixOperations::solveKramer( matrix, rightParts );
		}
		// если решений нет - прерываемся, иначе вычисляем начальные условия для следующей итерации
		if ( next.size() == 0 || ++iterCount > maxIterCount )
		{
			notFinded = true;
		}
		else
		{
			// вычисляем норму ( максимальную дельту )
			norm = calcNorm( next );
			// вычисляем начальные значения для следующей итерации
			for ( size_t i = 0; i < next.size(); i++ )
			{
				system.getUnknownVariables()[i].setValue( system.getUnknownVariables()[i].getValue() + next[i] );
			}
		}
	} while ( !notFinded && norm >= eps );

	NS_CORE InternalGearRatios ans;
	if ( !notFinded )
	{
		ans = geInternalGearRatioValuesFromSystem( system );
	}

	return ans;
}

int main(int argc, char** argv)
{
  int i, j, N, flag;
  Vector grid;
  Matrix b, e;
  double time, sum, h;
  int rank, size;
  int mpi_top_coords[2];
  int mpi_top_sizes[2];

  init_app(argc, argv, &rank, &size);

  N=atoi(argv[1]);

  // setup topology
  mpi_top_sizes[0] = mpi_top_sizes[1] = 0;
  MPI_Dims_create(size, 2, mpi_top_sizes);
  int periodic[2] = {0, 0};
  MPI_Comm comm;
  MPI_Cart_create(MPI_COMM_WORLD, 2, mpi_top_sizes, periodic, 0, &comm);
  MPI_Cart_coords(comm, rank, 2, mpi_top_coords);

  int* size1;
  int* displ1;
  int* size2;
  int* displ2;
  splitVector(N, mpi_top_sizes[0], &size1, &displ1);
  splitVector(N, mpi_top_sizes[1], &size2, &displ2);

  b = createMatrix(size1[mpi_top_coords[0]], size2[mpi_top_coords[1]]);
  for (j=0;j<b->cols;++j)
    for(i=0;i<b->rows;++i)
      b->data[j][i] = (j+displ2[mpi_top_coords[1]])*N+1+(i+displ1[mpi_top_coords[0]]);
  b->glob_rows = N;
  b->glob_cols = N;
  b->as_vec->comm = &comm;

  saveMatrixMPI(b, "meh.asc");

  freeMatrix(b);
  MPI_Comm_free(&comm);

  close_app();
  return 0;
}