int main () {
	struct Matrix *a = create_matrix(3, 3);
	struct Matrix *b;
	
	read_matrix(a);
	b = matrix_pow(a, 2);

	print_matrix(b);

	return 0;
}

int main(void) {
	int i, j;
	char s[MAXLEN], t[MAXLEN];		/* input strings */

	s[0] = t[0] = ' ';
	scanf("%s", &(s[1]));	scanf("%s", &(t[1]));

	printf("matching cost = %d \n", string_compare(s, t, strlen(s)-1,
	  strlen(t)-1));

	print_matrix(s, t, TRUE);
	printf("\n");
	print_matrix(s, t, FALSE);
	goal_cell(s, t, &i, &j);
	printf("%d %d\n", i, j);
	reconstruct_path(s, t, i, j);
	printf("\n");
	
	return 0;
}

int main() {
	
	Matrix *A = init_matrix(3, 3);
	Matrix *B = init_matrix(3, 2);
	
	(A->values)[0][0] = 0;
	(A->values)[0][1] = 1;
	(A->values)[0][2] = 2;
	(A->values)[1][0] = 3;
	(A->values)[1][1] = 4;
	(A->values)[1][2] = 5;
	(A->values)[2][0] = 1;
	(A->values)[2][1] = 1;
	(A->values)[2][2] = 1;
	
	(B->values)[0][0] = 0;
	(B->values)[0][1] = 1;
	(B->values)[1][0] = 2;
	(B->values)[1][1] = 3;
	(B->values)[2][0] = 4;
	(B->values)[2][1] = 5;
	
	print_matrix(A);
	print_matrix(B);
	
	Matrix *C = mmult(A, B);
	print_matrix(C);
	
			
	destroy_matrix(A);
	destroy_matrix(B);
	destroy_matrix(C);
	
	
	
	
	
	
	
	return 0;
}

void print_matrix_types() {

  double *random, *ones, *zeros, *identity, *tri;

  /* Allocate matrices */
  random = random_matrix(6, 3);
  identity = identity_matrix(5, 5);
  ones = ones_matrix(4, 2);
  zeros = zeros_matrix(2, 4);
  tri = lowerTri_matrix(5, 5);

  printf("\n\t\tMatrix Types\n");

  printf("6x3 Random Matrix\n");
  print_matrix(6, 3, random);
  printf("\n\n");

  printf("5x5 Identity Matrix\n");
  print_matrix(5, 5, identity);
  printf("\n\n");

  printf("4x2 Ones Matrix\n");
  print_matrix(4, 2, ones);
  printf("\n\n");

  printf("2x4 Zeros Matrix\n");
  print_matrix(2, 4, zeros);
  printf("\n\n");

  printf("5x5 Lower Triangular Matrix\n");
  print_matrix(5, 5, tri);
  printf("\n\n");

  /* deallocate memory */
  deallocate_matrix(random);
  deallocate_matrix(ones);
  deallocate_matrix(zeros);
  deallocate_matrix(identity);
  deallocate_matrix(tri);

}

int main(int argc, char* argv[]) {
    if (argc != 2) {
        printf("argc error!\n");
        return -1;
    }
    struct Matrix m;
    if (0 == init_matrix_from_file(&m, argv[1])) {
        print_matrix(&m);
        get_saddle(&m);
    }
    return 0;
}

int main(){
    double A[2][2] = {{4., 10.}, {1., 1.}};
    double inverseA[2][2];
    
    double determinant = A[1][1]*A[0][0] - A[1][0]*A[0][1];
    assert(determinant != 0);
    
    double inverse_det = pow(determinant, -1);
    
    inverseA[0][0] =  inverse_det * A[1][1];
    inverseA[1][0] = - inverse_det * A[1][0];
    inverseA[1][1] =  inverse_det * A[0][0];
    inverseA[0][1] = - inverse_det * A[0][1];
    
    
    std::cout << "Matrix A:\n";
    print_matrix(A);
    std::cout << "Inverse A:\n";
    print_matrix(inverseA);
    return 0;
}

void to_matrix(Circuit* c,gsl_matrix** A_ptr,gsl_vector** b_ptr)
{
	int w=c->ccount+c->vcount;
	(*A_ptr)=gsl_matrix_calloc(w,w);
	(*b_ptr)=gsl_vector_calloc(w);
	gsl_matrix * A=*A_ptr;
	gsl_vector * b=*b_ptr;
	if(A==NULL||b==NULL){
		printf("ERROR, UNABLE TO ALLOCATE MEMORY FOR MATRIX CONVERSION\n");
		if(A!=NULL){
			gsl_matrix_free(A);
			A=NULL;
		}
		if(b!=NULL){
			gsl_vector_free(b);
			b=NULL;
		}
		return;
	}
	int i;
	int j;
	for(i=0;i<c->ccount;i++)
	{
		Component* edge=c->components+i;
		for(j=0;j<c->vcount-1;j++){
			if(edge->A==c->vertices[j].id){
				gsl_matrix_set(A,j,i,-1.0);
			}
			if(edge->B==c->vertices[j].id){
				gsl_matrix_set(A,j,i,1.0);
			}
		}
	}
	gsl_matrix_view voltage_view=gsl_matrix_submatrix(A,c->vcount-1,c->ccount,c->ccount,c->vcount);
	gsl_matrix* voltage=&voltage_view.matrix;
	for(i=0;i<c->ccount;i++){
		for(j=0;j<c->vcount;j++){
			if(c->vertices[j].id==c->components[i].A)gsl_matrix_set(voltage,i,j,-1);
			if(c->vertices[j].id==c->components[i].B)gsl_matrix_set(voltage,i,j,1);
		}
		if(c->components[i].type==RESISTOR)
		{
			gsl_matrix_set(A,c->vcount+i-1,i,*(double *)c->components[i].data);
		}
		if(c->components[i].type==BATTERY)
		{
			gsl_vector_set(b,c->vcount+i-1,*(double *)c->components[i].data);
		}
	}
	gsl_matrix_set(A,w-1,w-1,1);
	print_matrix(voltage);
	print_vector(b);
}

int main() {
  
  igraph_t g;
  igraph_matrix_t m;
  
  igraph_small(&g, 0, IGRAPH_DIRECTED, 
	       0,1, 2,1, 2,0, 3,0, 
	       -1);
  
  igraph_matrix_init(&m, 0, 0);
  igraph_bibcoupling(&g, &m, igraph_vss_all());
  print_matrix(&m, stdout);
  
  igraph_cocitation(&g, &m, igraph_vss_all());
  print_matrix(&m, stdout);
  
  igraph_matrix_destroy(&m);
  igraph_destroy(&g);
  
  return 0;
}

int main(){
    float *a, *b, *c; // one dementional array
    int i, x, y;
    x = 3;
    y = 4;
    a = (float *)malloc(x*y*sizeof(float));
    for(i=0;i<x*y;++i){
        a[i] = i;
    }
    print_matrix(a,x,y);
    return 0;
}

void expDisk::load(host_state_t &hstate, const peyton::system::Config &cfg)
{
	// density distribution parameters
	l        = cfg.get("l");
	h        = cfg.get("h");

	assert(l > 0);
	assert(h > 0);
#if 0
	float z0 = cfg.get("z0");	// Solar offset from the Galactic plane
	Rg       = cfg.get("Rg");	// Distance to the Galactic center (assumed to be in l=0, b=0 direction)

	// compute the rotation of the Galactic plane, and cylindrical
	// coordinates of the Sun
	zsun = z0;
	rsun = sqrt(double(Rg*Rg) - double(z0*z0));
	asin = zsun / Rg;
	acos = rsun / Rg;
#else
	std::string ctr, orient;
	cfg.get(ctr,    "center",      "galplane");
	cfg.get(orient, "orientation", "galplane");
	load_transform(&T.x, M, ctr, orient, cfg);
#if 0
	std::cout << "translation = " << std::setprecision(10) << T.x << " " << T.y << " " << T.z << "\n";
	print_matrix(M);

	// Do some testing here...
	float3 v = { 0.f, 0.f, 0.f };
	v = transform(v, T, M);
	std::cout << std::setprecision(10) << v.x << " " << v.y << " " << v.z << "\n";
	abort();
#endif
#endif

	int userComp = cfg.get("comp");
	comp = componentMap.seqIdx(userComp);

	// set this to 0. for now, to allow LF normalization
	// even beyond the user-specified model cutoff
	r_cut2 = 0.f;

	// Load luminosity function
	std::string lffile;
	hstate.lf = load_lf(*this, cfg, lffile);

	// cutoff radius (default: no cutoff)
	cfg.get(r_cut2,  "rcut",   0.f);
	r_cut2 *= r_cut2;

	if(lffile.empty()) { lffile = "-"; }
	MLOG(verb1) << "Component " << userComp << " : " << "exponential disk {" << l << ", " << h << ", " << lffile << "}";
}

int main()
{
    Crosslist mat;
    creat_matrix(&mat);
    print_matrix(&mat);
    struct olnode *q;
    q = (&mat)->rhead[1];
    printf ("%d\n", q->e);
    q = q->right;
    printf ("%d\n", q->e);
    return 0;
}

/* Main program */
int main() {
        /* Locals */
        int m = M, n = N, lda = LDA, ldu = LDU, ldvt = LDVT, info, lwork;
        double wkopt;
        double* work;
        /* Local arrays */
        double s[N], u[LDU*M], vt[LDVT*N];
        double a[LDA*N] = {
            8.79,  6.11, -9.15,  9.57, -3.49,  9.84,
            9.93,  6.91, -7.93,  1.64,  4.02,  0.15,
            9.83,  5.04,  4.86,  8.83,  9.80, -8.99,
            5.45, -0.27,  4.85,  0.74, 10.00, -6.02,
            3.16,  7.98,  3.01,  5.80,  4.27, -5.31
        };
        /* Executable statements */
        printf( " DGESVD Example Program Results\n" );
        /* Query and allocate the optimal workspace */
        lwork = -1;
        dgesvd( "All", "All", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt, &lwork,
         &info );
        lwork = (int)wkopt;
        work = (double*)malloc( lwork*sizeof(double) );
        /* Compute SVD */
        dgesvd( "All", "All", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork,
         &info );
        /* Check for convergence */
        if( info > 0 ) {
                printf( "The algorithm computing SVD failed to converge.\n" );
                exit( 1 );
        }
        /* Print singular values */
        print_matrix( "Singular values", 1, n, s, 1 );
        /* Print left singular vectors */
        print_matrix( "Left singular vectors (stored columnwise)", m, n, u, ldu );
        /* Print right singular vectors */
        print_matrix( "Right singular vectors (stored rowwise)", n, n, vt, ldvt );
        /* Free workspace */
        free( (void*)work );
        exit( 0 );
} /* End of DGESVD Example */

int main() {
    double A[16] = { -1, 0, 1, 8,
                     1, 0, 1, 4,
                     -1, 0, 1, 2,
                     1, 1, 1, 1};
    double R[16];
    
    matrix_invert(4, A, (double *)R);

    print_matrix(4, 4, R);

    return 0;
}

int main(int argc,char *argv[])
{
if(argc!=2)
{
printf("No of rows not specified \n");
return 1;
}
struct timeval tv,tv1;
nrows=atoi(argv[1]);
ncols=nrows;
mapsize=nrows*ncols*sizeof(float);
fill_matrix("./matrixA");
fill_matrix("./matrixB");
fill_matrix("./matrixC");
printf("Matrix filled \n");
ma=map_matrix("./matrixA",MAP_RDONLY,nrows,ncols,&fda);
mb=map_matrix("./matrixB",MAP_RDONLY,nrows,ncols,&fdb);
mc=map_matrix("./matrixC",MAP_RDWR,nrows,ncols,&fdc);
gettimeofday(&tv,NULL);
par_vector_vector_multiply();
gettimeofday(&tv1,NULL);
double tm=tv1.tv_sec-tv.tv_sec+(tv1.tv_usec-tv.tv_usec)*0.000001;
printf("%lf , %lf \n",tv1.tv_sec+0.0,tv.tv_sec+0.0);
printf("\n time taken is %lf  seconds \n",tm);

printf("Matrix A is ..\n");
print_matrix(ma);
printf("Matrix B is ..\n");
print_matrix(mb);
printf("Matrix C is ..\n");
print_matrix(mc);
munmap(ma,mapsize);
munmap(mb,mapsize);
munmap(mc,mapsize);
close(fda);
close(fdb);
close(fdc);
return 0;
}

int main(void)
{
	int i, j, k;
	float a[m][n];
	for(i=0; i<m; i++)
	{
		for(j=0; j<n; j++)
		{
			a[m][n] = cos(((float)n/(m+1)) * M_PI);
		}
	}
	print_matrix(a[m][m], m, n);
}

int main(int argc, char** argv)

{
    int D=10;
    srand(time(0));
    double *A;
    double *B;
    init_matrix(&A,D,D);
    init_matrix(&B,D,1);
    printf("The first input matrix is:\n");
    print_matrix(A,D,D);
    printf("The Second input matrix is:\n");
    print_matrix(B,D,1);
    double *x = malloc(D*sizeof(double));
    solve_system(A, D, B, x);
    printf("The result is:\n");
    print_matrix(x,D,1);
    printf("\n");


    return 0;
}

int
abg_to_xyz(int argc, char *argv[])
{
  PBASIS p;
  XVBASIS xv;
  double  **covar_abg, **covar_xyz, **derivs;

  int	i,j;

/*  if (argc!=1) print_help(abg_to_xyz_help); */

  covar_abg = dmatrix(1,6,1,6);
  covar_xyz = dmatrix(1,6,1,6);
  derivs = dmatrix(1,6,1,6);

  /* echo the command line to output */
  printf("#");
  for (i=0; i<argc; i++) printf(" %s",argv[i]);
  {
#include <time.h>
    time_t timettt;
    time(&timettt);
    /* note that ctime returns string with newline at end */
    printf("\n#---%s",ctime(&timettt));
  }

  if (read_abg(NULL,&p,covar_abg)) {
    fprintf(stderr, "Error in input alpha/beta/gamma file\n");
    exit(1);
  }

  /* Transform the orbit basis and get the deriv. matrix */
  pbasis_to_bary(&p, &xv, derivs);

  /* Map the covariance matrix to new basis */
  covar_map(covar_abg, derivs, covar_xyz,6,6);

  /* Print out the results, with comments */
  printf("#     x              y             z           "
	 "vx          vy            vz\n");
  printf("%12.8f %12.8f %12.8f %12.8f %12.8f %12.8f\n",
	 xv.x, xv.y, xv.z, xv.xdot, xv.ydot, xv.zdot);
  printf("# covariance matrix:\n");
  print_matrix(stdout,covar_xyz,6,6);
  printf("# Position is for JD:\n%.5f\n",jd0);

  free_dmatrix(covar_abg,1,6,1,6);
  free_dmatrix(covar_xyz,1,6,1,6);
  free_dmatrix(derivs,1,6,1,6);
  exit(0);
}

void prob2a(){
  std::cout << "Problem 2 Part A ======== " << std::endl;

  Matrix points1(10000,std::vector<double>(2));
  Matrix points2(10000,std::vector<double>(2));
  loadFile(points1, "sample_data/output1");
  loadFile(points2, "sample_data/output3");

  std::vector<double> mean1 = getSampleMean(points1);
  std::vector<double> mean2 = getSampleMean(points2);

  std::cout << "sample_mean1 = ";
  print_vec(mean1);
  std::cout << "sample_mean2 = ";
  print_vec(mean2);

  Matrix cov1 = getSampleVar(points1, mean1);
  Matrix cov2 = getSampleVar(points2, mean2);

  std::cout << "sample_cov1 = ";
  print_matrix(cov1);
  std::cout << "sample_cov2 = ";
  print_matrix(cov2);

  std::cout << "With estimated params: " << std::endl;
  QuadraticDiscriminant classifier1(mean1, mean2, cov1, cov2);
  calcError(classifier1,points1,points2, "sample_data/labels1");

  std::cout << "Given params: " << std::endl;
  mean1 = {1.0,1.0};
  mean2 = {6.0,6.0};
  cov1 = {{2.0,0},{0,2.0}};
  cov2 = {{4.0,0},{0,8.0}};

  QuadraticDiscriminant classifier2(mean1, mean2, cov1, cov2);
  calcError(classifier2,points1,points2, "sample_data/labels1");
  std::cout << std::endl;
}

/* {{{ check */
int check(matrix *m) {
    int i, j;
    REAL temp;
    init_empty_matrix(&lower, m->rows, m->cols);
    init_empty_matrix(&upper, m->rows, m->cols);
    for (i=0; i<m->rows; i++) {
        for (j=0; j<m->cols; j++) {
            if (i>j) {
                temp = GET(m,i,j);
                SET(&lower,i,j,temp);
            }
            if (i<=j) {
                temp = GET(m,i,j);
                SET(&upper,i,j,temp);
            }
        }
        SET(&lower,i,i,(REAL)1.0L);
    }
    mult_matrices(&mult, &lower, &upper);
    if (dbg) {
        printf("lower:\n");
        print_matrix(&lower);
        printf("upper:\n");
        print_matrix(&upper);
        printf("multiplied:\n");
        print_matrix(&mult);
    }
    for (i=0; i<m->rows; i++) {
        for (j=0; j<m->cols; j++) {
            if (!EQUAL(GET(&mult,i,j),GET(&ref,i,j))) {
                printf("(%d,%d): " FMT " != " FMT "\n",
                        i, j, GET(&mult,i,j), GET(&ref,i,j));
                return EXIT_FAILURE;
            }
        }
    }
    return EXIT_SUCCESS;
} /* }}} */

int inv_test(u8 * in, u8 * inv, u8 * sav, int n)
{
	memcpy(sav, in, n * n);

	if (gf_invert_matrix(in, inv, n)) {
		printf("Given singular matrix\n");
		print_matrix(sav, n);
		return -1;
	}

	matrix_mult(inv, sav, in, n);

	if (is_ident(in, n)) {
		printf("fail\n");
		print_matrix(sav, n);
		print_matrix(inv, n);
		print_matrix(in, n);
		return -1;
	}
	putchar('.');

	return 0;
}