/* matroot.c - calculate the square root of a matrix
 *
 * Copyright (C) 2004  Jochen Voss.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * $Id: matroot.c 7208 2006-09-03 16:11:08Z voss $
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#ifdef USECBLAS
#include <cblas.h>
#endif

int N = 7;


static void
set_entry(double *A, int i, int j, double val)
{
  A[j*N+i] = val;
}

static double
get_entry(const double *A, int i, int j)
{
  return A[j*N+i];
}

void  my_dgemm(int size, double *a, double *b, double *c) 
{
  int i, j, k;
    for (i=0; i<size; i++)
      for (j=0; j < size; j++)
        for (k=0; k<size; k++) {
          *(c+i*size + j) = *(c+i*size + j) + *(a+i*size + k) *  *(b +k*size + j);
	}
}


int
main( int argc, char *argv[])
{
  double *A, *B, *C, *W, *Z, *WORK;
  int *ISUPPZ, *IWORK;
  int  i, j;
  int  M;

  if (argc != 3) {
    printf("Usage: a.out matrix_size iter\n");
    exit(0);
  }

  N= atoi(argv[1]);
  M = atoi(argv[2]);
  /* allocate and initialise the matrix */
  A = malloc(N*N*sizeof(double));
  for (i=0; i<N; ++i) {
    for (j=0; j<i-1; ++j) {
      set_entry(A, i, j, 0);
    }
  }
  for (i=0; i<N-1; ++i)  set_entry(A, i+1, i, -1);
  for (i=1; i<N-1; ++i)  set_entry(A, i, i, 2);
  set_entry(A, 0, 0, 1);
  set_entry(A, N-1, N-1, 1);

  /* allocate space for the output parameters and workspace arrays */
  W = malloc(N*sizeof(double));
  Z = malloc(N*N*sizeof(double));
  ISUPPZ = malloc(2*N*sizeof(int));
  WORK = malloc(26*N*sizeof(double));
  IWORK = malloc(10*N*sizeof(int));

  /* allocate and initialise a new matrix B=Z*D */
  B = malloc(N*N*sizeof(double));
  for (j=0; j<N; ++j) {
    double  lambda=sqrt(W[j]);
    for (i=0; i<N; ++i) {
      set_entry(B, i, j, get_entry(Z,i,j)*lambda);
    }
  }

  C = malloc(N*N*sizeof(double));

  for (i=0; i<M; i++) {

#ifdef USECBLAS 
  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, N, N, N,
              1, B, N, A, N, 0, C, N);
#else 
  my_dgemm(N, B, A, C);
#endif
  }

  return 0;
}