static void
compute_lhs(const struct efp *efp, double *c, bool conj)
{
	size_t n = 3 * efp->n_polarizable_pts;

	for (size_t i = 0, offset_i = 0; i < efp->n_frag; i++) {
	for (size_t ii = 0; ii < efp->frags[i].n_polarizable_pts; ii++, offset_i++) {
	for (size_t j = 0, offset_j = 0; j < efp->n_frag; j++) {
	for (size_t jj = 0; jj < efp->frags[j].n_polarizable_pts; jj++, offset_j++) {
		if (i == j) {
			if (ii == jj)
				copy_matrix(c, n, offset_i, offset_j, &mat_identity);
			else
				copy_matrix(c, n, offset_i, offset_j, &mat_zero);

			continue;
		}

		const struct polarizable_pt *pt_i = efp->frags[i].polarizable_pts + ii;
		mat_t m = get_int_mat(efp, i, j, ii, jj);

		if (conj)
			m = mat_trans_mat(&pt_i->tensor, &m);
		else
			m = mat_mat(&pt_i->tensor, &m);

		mat_negate(&m);
		copy_matrix(c, n, offset_i, offset_j, &m);
	}}}}
}

    void
    factor_inner (const Ordinal m,
		  const Ordinal n,
		  Scalar R[],
		  const Ordinal ldr,
		  Scalar A[],
		  const Ordinal lda,
		  Scalar tau[],
		  Scalar work[])
    {
      const Ordinal numRows = m + n;

      A_buf_.reshape (numRows, n);
      A_buf_.fill (Scalar(0));
      // R might be a view of the upper triangle of a cache block, but
      // we only want to include the upper triangle in the
      // factorization.  Thus, only copy the upper triangle of R into
      // the appropriate place in the buffer.
      copy_upper_triangle (n, n, &A_buf_(0, 0), A_buf_.lda(), R, ldr);
      copy_matrix (m, n, &A_buf_(n, 0), A_buf_.lda(), A, lda);

      int info = 0;
      lapack_.GEQR2 (numRows, n, A_buf_.get(), A_buf_.lda(), tau, work, &info);
      if (info != 0)
	throw std::logic_error ("TSQR::CombineDefault: GEQR2 failed");

      // Copy back the results.  R might be a view of the upper
      // triangle of a cache block, so only copy into the upper
      // triangle of R.
      copy_upper_triangle (n, n, R, ldr, &A_buf_(0, 0), A_buf_.lda());
      copy_matrix (m, n, A, lda, &A_buf_(n, 0), A_buf_.lda());
    }

/***********************************************************************
 * Copy a motif from one place to another.
 ***********************************************************************/
void copy_motif
  (MOTIF_T* source,
   MOTIF_T* dest)
{
  strcpy(dest->id, source->id);
  strcpy(dest->id2, source->id2);
  dest->length = source->length;
  dest->alph_size = source->alph_size;
  dest->ambigs = source->ambigs;
  dest->evalue = source->evalue;
  dest->num_sites = source->num_sites;
  dest->complexity = source->complexity;
  dest->trim_left = source->trim_left;
  dest->trim_right = source->trim_right;
  if (source->freqs) {
    // Allocate memory for the matrix.
    dest->freqs = allocate_matrix(dest->length, dest->alph_size + dest->ambigs);
    // Copy the matrix.
    copy_matrix(source->freqs, dest->freqs);
  } else {
    dest->freqs = NULL;
  }
  if (source->scores) {
    // Allocate memory for the matrix. Note that scores don't contain ambigs.
    dest->scores = allocate_matrix(get_num_rows(source->scores), get_num_cols(source->scores));
    // Copy the matrix.
    copy_matrix(source->scores, dest->scores);
  } else {
    dest->scores = NULL;
  }
  copy_string(&(dest->url), source->url);
}

/** Matrix inversion via the Gauss-Jordan algorithm. */
static elem_t* invert_matrix(const elem_t* const a, const int n)
{
  int i, j;
  elem_t* const inv = new_matrix(n, n);
  elem_t* const tmp = new_matrix(n, 2*n);
  copy_matrix(a, n, n, 0, n, 0, n, tmp, n, 2 * n, 0, n, 0, n);
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      tmp[i * 2 * n + n + j] = (i == j);
  gj(tmp, n, 2*n);
  copy_matrix(tmp, n, 2*n, 0, n, n, 2*n, inv, n, n, 0, n, 0, n);
  delete_matrix(tmp);
  return inv;
}

matrix_t *
symm_pivot_ge_matrix (matrix_t * a, int *row_per)
{
  matrix_t *c = copy_matrix (a);
  if (c != NULL) {
    int i, j, k;
    int cn = c->cn;
    int rn = c->rn;
    double *e = c->e;
    for (i = 0; i < rn; i++)
      row_per[i] = i;
    for (k = 0; k < rn - 1; k++) {      /* eliminujemy (zerujemy) kolumnę nr k */
      int piv = k;              /* wybór eleemntu dominującego - maks. z k-tej kol., poniżej diag */
      for (i = k + 1; i < rn; i++)
        if (fabs (*(e + i * cn + k)) > fabs (*(e + piv * cn + k)))
          piv = i;
      if (piv != k) {           /* jeśli diag. nie jest pivtem - wymień wiersze */
        int tmp;
        xchg_rows (c, piv, k);
        xchg_cols (c, piv, k);
        tmp = row_per[k];
        row_per[k] = row_per[piv];
        row_per[piv] = tmp;
      }
      for (i = k + 1; i < rn; i++) {    /* pętla po kolejnych
                                           wierszach poniżej diagonalii k,k */
        double d = *(e + i * cn + k) / *(e + k * cn + k);
        for (j = k; j < cn; j++)
          *(e + i * cn + j) -= d * *(e + k * cn + j);
      }
    }
  }
  return c;
}

MATRIX *invert_ltriangle_matrix(MATRIX *L)
/* Given a lower triangular matrix L, computes inverse L^-1 */
{
  int i,j,k,n;
  double sum;
  MATRIX *I;
  
  if(L->m != L->n) {
    printf("ERROR: Matrix not quadratic. Cannot invert triangular matrix!\n");
    exit(1);
  }
  n=L->n;
  I=copy_matrix(L);

  for (i=0;i<n;i++) {
    I->element[i][i]=1.0/L->element[i][i];
    for (j=i+1;j<n;j++) {
      sum=0.0;
      for (k=i;k<j;k++) sum -= I->element[j][k]*I->element[k][i];
      I->element[j][i]=sum/L->element[j][j];
    }
  }

  return(I);
}

/************************************************************************
 * See .h file for description.
 ************************************************************************/
void copy_mhmm
  (MHMM_T*        an_mhmm,
   MHMM_T*        new_mhmm)
{
  int i_state;

  /* Copy the top-level data. */
  new_mhmm->type = an_mhmm->type;
  new_mhmm->log_odds = an_mhmm->log_odds;
  new_mhmm->num_motifs = an_mhmm->num_motifs;
  new_mhmm->num_states = an_mhmm->num_states;
  new_mhmm->num_spacers = an_mhmm->num_spacers;
  new_mhmm->spacer_states = an_mhmm->spacer_states;
  new_mhmm->alph = alph_hold(an_mhmm->alph);
  new_mhmm->background = allocate_array(alph_size_full(an_mhmm->alph));
  copy_array(an_mhmm->background, new_mhmm->background);
  copy_string(&(new_mhmm->description), an_mhmm->description);
  copy_string(&(new_mhmm->motif_file), an_mhmm->motif_file);
  copy_string(&(new_mhmm->sequence_file), an_mhmm->sequence_file);
  // FIXME: Copy hot states array.
  new_mhmm->num_hot_states = an_mhmm->num_hot_states;

  /* Copy each state. */
  for (i_state = 0; i_state < an_mhmm->num_states; i_state++) {
    copy_state(&(an_mhmm->states[i_state]),
               &(new_mhmm->states[i_state]));
  }

  /* Copy the transition matrix. */
  copy_matrix(an_mhmm->trans, new_mhmm->trans);
}

/******************************************************************************
  * This function returns a copy of the matrix for time t of a transition 
  * matrix array. Returns NUL is no matrix is found matching time t.
  * Caller is responsible for free the returned copy.
******************************************************************************/
MATRIX_T* get_substmatrix_for_time(
  SUBSTMATRIX_ARRAY_T *a, 
  double t
) {

  assert(a != NULL);
  // Create a copy of our matrix
  int i;
  MATRIX_T* src = NULL;
  for (i = 0; i < a->length; i++) {
    if (a->t[i] == t) {
      src = a->substmatrix[i];
      break;
    }
  }
  MATRIX_T* dst = NULL;
  if (src != NULL) {
    int num_rows = get_num_rows(src);
    int num_cols = get_num_cols(src);
    dst = allocate_matrix(num_rows, num_cols);
    copy_matrix(src, dst);
  }

  return dst;
}

matrix * WEKAgetSampleAccuracy(struct hash * config)
{
int fold, folds = foldsCountFromDataDir(config);
int split, splits = splitsCountFromDataDir(config);

matrix * accuracies = NULL;
for(split = 1; split <= splits; split++)
	{
	for(fold = 1; fold <= folds; fold++)
	    {
		matrix * accuraciesInFold = WEKApopulateAccuracyMatrix(config,split, fold);

	    //add the accuracies to the running totals
	    if(split == 1 && fold == 1) 
	        accuracies = copy_matrix(accuraciesInFold);
	    else
	        add_matrices_by_colLabel(accuracies, accuraciesInFold);
	
	    //clean up
	    free_matrix(accuraciesInFold);
		}
	}
//normalize accuracies over number of splits and folds
int i;
for(i = 0; i < accuracies->cols; i++)
    {
    if(accuracies->graph[0][i] != NULL_FLAG)
        accuracies->graph[0][i] = (accuracies->graph[0][i] / ((folds-1) * splits));
    if(accuracies->graph[1][i] != NULL_FLAG)
        accuracies->graph[1][i] = (accuracies->graph[1][i] / (1 * splits));
    }
return accuracies;
}

/* should output:
 * print_matrix:
 *  3  1 -1
 *  0 -2  0
 *  5  0  0
 * matrix minor:
 *  0  0
 *  5  0
 * ineff_det: -10
 * eff_det: 10 
 * lu decomposition is not unique, output can be checked manually
 * invert_lower_tri_matrix: 
 * 1.000000 0.000000 0.000000 
 * 0.000000 1.000000 0.000000 
 * -1.666667 -0.833333 1.000000 
 * invert_upper_tri_matrix: 
 * 0.333333 0.166667 0.200000 
 * 0.000000 -0.500000 0.000000 
 * 0.000000 -0.000000 0.600000 
 * invert_matrix: 
 * 0.000000 0.000000 0.200000 
 * 0.000000 -0.500000 0.000000 
 * -1.000000 -0.500000 0.600000 */
void test_matrix_functions() {
	matrix *a = identity_matrix(3);
	a->entries[0][0] = 3.0;
	a->entries[0][1] = 1.0;
	a->entries[0][2] = -1.0;
	a->entries[1][1] = -2.0;
	a->entries[2][0] = 5.0;
	a->entries[2][2] = 0.0;
	printf("print_matrix: \n");
	print_matrix(*a);
	printf("matrix_minor: \n");
	print_matrix(*matrix_minor(*a, 1));
	printf("ineff_det: %lf\n", (double) ineff_det(*a));
	printf("eff_det: %lf\n", (double) eff_det(*a));
	matrix *a_cpy = copy_matrix(*a);
	printf("lu_decomp: \n");
	matrix **pl = lu_decomp(a_cpy, (int*) NULL);
	print_matrix(*pl[0]);	// prints p
	print_matrix(*pl[1]);	// prints l
	print_matrix(*a_cpy);	// prints u
	printf("invert_lower_tri_matrix: \n");
	print_matrix(*invert_lower_tri_matrix(*pl[1]));
	printf("invert_upper_tri_matrix: \n");
	print_matrix(*invert_upper_tri_matrix(*a_cpy));
	printf("invert_matrix: \n");
	print_matrix(*invert_matrix(*a));
}

void multiply_top_stack(matrix * m){
	matrix *copy = new matrix();
	matrix *ptr = peek_stack();
	copy_matrix(ptr,copy);	
	matrix_multiply(copy,m,ptr);
	delete(copy);
}

/*-------------- void matrix_mult() --------------
Inputs:  struct matrix *a
         struct matrix *b 
Returns: 

a*b -> b
*/
void matrix_mult(struct matrix *a, struct matrix *b) {

  int r;
  double t;
  int i, j;

  if (a->cols != b->rows){
    printf ("Improperly sized matrices!\n");
    return;
  }

  struct matrix *c = new_matrix(b->rows, b->cols);
   copy_matrix(b, c);
   c->lastcol = b->lastcol;

  *b = *new_matrix(a->rows, b->cols);
  for (i = 0; i < a->rows; i++){
    for(j = 0; j < c->cols; j++){
      r = 0;
      t = 0;
      while(r < a->rows){
	t+= (a-> m[i][r])*(c-> m[r][j]);
      	r++;
      }
      b->m[i][j] = t;
    }
  }
  free_matrix(c);
}

void push(stack *hay){
  hay -> top ++;
  hay -> matrices[hay -> top] = *new_matrix(4, 4);
  copy_matrix(&hay -> matrices[hay -> top - 1], &hay -> matrices[hay -> top]);
  //printf("pushed: %d\n", hay -> top);
  //print_matrix(&hay -> matrices[hay -> top]);
}

 /// \brief Copy constructor.
 ///
 /// We need an explicit copy constructor, because otherwise the
 /// default copy constructor would override the generic matrix
 /// view "copy constructor" below.
 Matrix (const Matrix& in) :
     nrows_ (in.nrows()),
     ncols_ (in.ncols()),
     A_ (verified_alloc_size (in.nrows(), in.ncols()))
 {
     if (! in.empty())
         copy_matrix (nrows(), ncols(), get(), lda(), in.get(), in.lda());
 }

static bool set_attribute_matrix(const Transform &tfm, TypeDesc type, void *val)
{
  if (type == TypeDesc::TypeMatrix) {
    copy_matrix(*(OSL::Matrix44 *)val, tfm);
    return true;
  }

  return false;
}

extern MATRIX_T* gen_pam_matrix(
  ALPH_T alph,                  /* alphabet */
  int dist,			/* PAM distance */
  BOOLEAN_T logodds		/* true: generate log-odds matrix 
				   false: generate target frequency matrix 
				*/
)
{
  assert(alph == DNA_ALPH || alph == PROTEIN_ALPH);
  int i, j;
  MATRIX_T *matrix, *mul;
  BOOLEAN_T dna = (alph == DNA_ALPH);
  double *pfreq = dna ? pam_dna_freq : pam_prot_freq;	// standard frequencies
  int alen = alph_size(alph, ALPH_SIZE);  // length of standard alphabet
  double factor = dist < 170 ? 2/log(2) : 3/log(2);	// same as in "pam" Version 1.0.6

  /* create the array for the joint probability matrix */
  matrix = allocate_matrix(alen, alen);
  mul = allocate_matrix(alen, alen);

  /* initialize the matrix: PAM 1:
     due to roundoff, take the average of the two estimates of the joint frequency
     of i and j as the joint, then compute the conditionals for the matrix
  */
  for (i=0; i<alen; i++) {
    for (j=0; j<=i; j++) {
      double vij = dna ? trans[i][j] : dayhoff[i][j];
      double vji = dna ? trans[j][i] : dayhoff[j][i];
      double joint = ((vij * pfreq[j]) + (vji * pfreq[i]))/20000;/* use average to fix rndoff */
      set_matrix_cell(i, j, joint/pfreq[j], matrix);
      if (i!=j) set_matrix_cell(j, i, joint/pfreq[i], matrix);
    }
  }

  /* take PAM matrix to desired power to scale it */ 
  copy_matrix(matrix, mul);
  for (i=dist; i>1; i--) {
    MATRIX_T *product = matrix_multiply(matrix, mul);
    SWAP(MATRIX_T*, product, matrix)
    free_matrix(product);
  } 
  free_matrix(mul);

  /* convert to joint or logodds matrix:
     target:  J_ij = Pr(i,j) = Mij pr(j) 
     logodds: L_ij = log (Pr(i,j)/(Pr(i)Pr(j)) = log (Mij Pr(j)/Pr(i)Pr(j)) = log(Mij/pr(i)) 
  */
  for (i=0; i<alen; i++) {
    for (j=0; j<alen; j++) {
      double vij = get_matrix_cell(i, j, matrix);
      vij = logodds ? nint(factor * log((vij+EPSILON)/pfreq[i])) : vij * pfreq[j];
      set_matrix_cell(i, j, vij, matrix);
    }
  }

  return matrix;
} /* gen_pam_matrix */

void			transform_matrix(double **matrix1, double **matrix2)
{
	double		**cpy;

	cpy = malloc_matrix(4, 4);
	matrix_multiply(matrix1, matrix2, cpy);
	copy_matrix(matrix1, cpy);
	free_matrix(cpy, 4);
}

int main(int argc, char* argv[])
{
    if(argc!=5){
	fprintf(stderr,"USAGE:\n\t%s <input file> <num steps> <output file> <max number of threads> \n", argv[0]);
	fprintf(stderr, "EXAMPLE:\n\t%s random.txt 500 result.txt 8\n", argv[0]);
	exit(EXIT_FAILURE);
    }

    char* init_fname = argv[1];
    sscanf(argv[2],"%d", &steps);
    char* out_fname = argv[3];
    int max_n_threads = atoi(argv[4]);

    long time_vs_threads[max_n_threads+1];
    int tmp[max_n_threads];
    threads_arg=tmp;

    char** init_matrix = dump_matrix_file(init_fname);
    char** zero_matrix = allocate_matrix(width, height);
    matrix_from = allocate_matrix(width, height);
    matrix_to = allocate_matrix(width, height);

    // run by different number of threads, and measure the time
    for(n_threads=1; n_threads <= max_n_threads; n_threads++){
	copy_matrix(init_matrix, matrix_from, width+2, height+2);
	copy_matrix(zero_matrix, matrix_to, width+2, height+2);
	
	time_vs_threads[n_threads] = walltime_of_threads(n_threads);
	printf("walltime when run by %d threads: %d microseconds\n", n_threads, time_vs_threads[n_threads]);
    }

    save_matrix_file(out_fname, matrix_from,width,height);

    char time_fname[20];
    sprintf(time_fname,"time-%d-%d.txt",width,height);
    save_vector(time_fname, time_vs_threads+1, max_n_threads);
   
    free_matrix(zero_matrix,width,height);
    free_matrix(init_matrix,width,height);
    free_matrix(matrix_from,width,height);
    free_matrix(matrix_to,width,height);
    return 0;
}

int matrix_invert__inline(double **initial, int rows, int cols) {
	double **result;
	int retval;

	new_matrix(&result, rows, rows);
	retval = matrix_invert(initial, rows, rows, result);
	copy_matrix(result, rows, rows, initial);
	free_matrix(&result, rows, rows);
	return retval;
}

static int my_reverse (double **a, double **b, unsigned int n)
{
   copy_matrix (b, a, n);
   if (-1 == JDM_ludecomp_inverse (b, n))
     {
	return -1;
     }
   
   return 1;
}