#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <lapacke.h>


// BLAS functions prototypes:

void drotg_(double *DA, double *DB, double *C, double *S);

void drot_(const int *N, double*DX, const int *INCX, double *DY, 
	const int *INCY, const double *C, const double *S);

double dasum_(const int *n, const double *X, const int *incx);

double dnrm2_(const int *n, const double *X, const int *incx);

/*
   LAPACKE functions prototypes:
   
   lapack_int LAPACKE_dgetrf 	( 	int  	matrix_order,
		lapack_int  	m,
		lapack_int  	n,
		double *  	a,
		lapack_int  	lda,
		lapack_int *  	ipiv 
	) 	
	
	lapack_int LAPACKE_dgecon 	( 	int  	matrix_order,
		char  	norm,
		lapack_int  	n,
		const double *  	a,
		lapack_int  	lda,
		double  	anorm,
		double *  	rcond 
	) 	
	
   double LAPACKE_dlange 	( 	int  	matrix_order,
		char  	norm,
		lapack_int  	m,
		lapack_int  	n,
		const double *  	a,
		lapack_int  	lda 
	) 	
	
	
	
	*/

#define index2D(i,j,LD1) i + (j*LD1)

double MAXVAL(int s, double *a)
{
   double v;
   int e;

   v = a[0];

   for ( e = 0; e < s; e++ )
   {
       if ( v < a[e] ) v = a[e];
   
   }
   
   return(v);
}
 
double MatrixNorm(const char *Norm, const int *N, const double *A, 
	const int *LDA, int *INFO) 
{
   double *TV, *ta, anorm;
   int i, j, rc, len, one;
   
   *INFO = 0;
   anorm = 0.0;
   if ( Norm[0] == '1' ) {
   	   // Column wise
		// Maximum column sum
			/// Code here
	} else if ( Norm[0] == '2' ) {
		/// square root of sum of squares
	} else if ( Norm[0] == 'i' || Norm[0] == 'I' ) {
	  // Row wise
	     // Maximum row sum
                /// Code here
	   
	} else {
    	*INFO = -1;
    }
	return(anorm);
}

double random_number() {
	/* random numbers 0 <= n < 1 */
	double r;
	
	r = (double)rand() / (double)RAND_MAX;
	
	return(r);
}


int main(  int argc, char *argv[] )
{ 
//
// Computes the Condition Number of a matrix A[LA,N] 
//
   double *A;    // A[LA,N] matrix
   int *Ipiv;  // workspace
   double x, y, xl, xr, yl, yu, dx, dy, Eps, Solution;
   double a1norm, a2norm, aInorm, da1norm, da2norm, daInorm;
   double r1cond, rIcond, rmin, rmax, v;
   int N, Np1, Nm1, LA, np, pos;
   int i, j, rc;
   lapack_int info;
   clock_t c0, c1;
   
   // Input reading with format:
   // LA N
   // A(1,1)    A(1,2)   A(1,3)   ... A(1,N)
   //		.	.	.
   // A(LA,1)  A(LA,2) A(LA,3) ... A(LA,N)
   //
   printf(" LA, N ? ");   
   fscanf(stdin,"%d, %d",&LA,&N);

   if ( LA < 0 || LA > 1000 ) LA = 1000;
   if ( N < 0 || N > 1000 ) N = 1000;
   printf(" LA = %d, N = %d\n",LA,N);
//   
   A = malloc( sizeof((double)1.0)*LA*N);
   if ( A == NULL ) {
   	  fprintf(stderr,"Error allocating A matrix\n");
      return(-1);
   }
//    

   Ipiv = malloc( sizeof((int)1)*LA);
   if ( Ipiv == NULL ) {
   	  fprintf(stderr,"Error allocating Ipiv\n");
      return(-1);
   } 

//
//  Fill A matrix 
//
   for ( i = 0; i < LA; i++ ) {
      for ( j = 0; j < N; j++ ) {
      	 if ( i == j) 
      	 	 A[index2D(j,i,LA)] = (double)1.0;
      	 else
      	 	 A[index2D(j,i,LA)] = random_number()/(double)1000.0;
   	  }
   }

   if ( LA <= 5 ) {
      printf(" A = \n");
      for ( i = 0; i < LA; i++ ) {
         for ( j = 0; j < N; j++ ) {
   	        printf("%lf ",A[index2D(j,i,LA)]);
      	  }
      	  printf("\n");
      }
   }

   printf("Computing matrix norm...\n"); 
   c0 = clock();
   a1norm = MatrixNorm("1",&LA,A,&LA,&rc);
   if ( rc != 0 ) {
   	  fprintf(stderr," After MatrixNorm('1') INFO = %d\n",rc);
   }
   a2norm = MatrixNorm("2",&LA,A,&LA,&rc);
   if ( rc != 0 ) {
   	  fprintf(stderr," After MatrixNorm('2') INFO = %d\n",rc);
   }
   aInorm = MatrixNorm("I",&LA,A,&LA,&rc);
   if ( rc != 0 ) {
   	  fprintf(stderr," After MatrixNorm('I') INFO = %d\n",rc);
   }
   c1 = clock();
   printf("Seconds for computing matrix norms: %lf\n",(double)(c1-c0)/(double)CLOCKS_PER_SEC);
   
   printf("Computing matrix norm with DLANGE ...\n"); 
   c0 = clock();
   da1norm = /// call LAPACKE_dlange properly (matrix_order=LAPACK_COL_MAJOR)
   da2norm = /// call LAPACKE_dlange properly (matrix_order=LAPACK_COL_MAJOR)
   daInorm = /// call LAPACKE_dlange properly (matrix_order=LAPACK_COL_MAJOR)
   c1 = clock();
   printf("Seconds for computing matrix norms with LAPACK: %lf\n",(double)(c1-c0)/(double)CLOCKS_PER_SEC);
   
   printf("Computing matrix factorization...\n"); 
   c0 = clock();
   for ( i = 0; i < LA; i++ ) Ipiv[i] = 0.0;
   
   /// call LAPACKE_dgetrf properly (matrix_order=LAPACK_COL_MAJOR)
   if ( info != 0 ) {
   	  fprintf(stderr," After DGETRF INFO = %d\n",info);
   }
   c1 = clock();
   printf("Seconds for matrix factorization: %lf\n",(double)(c1-c0)/(double)CLOCKS_PER_SEC);
  
   
   if ( LA <= 5 ) {
      printf(" Factorized A = \n");
      for ( i = 0; i < LA; i++ ) {
         for ( j = 0; j < N; j++ ) {
   	        printf("%lf ",A[index2D(j,i,LA)]);
      	  }
      	  printf("\n");
      }
   }
   
   printf("Computing Condition number:\n");  
   c0 = clock();
   /// call LAPACKE_dgecon with norm 1 properly (matrix_order=LAPACK_COL_MAJOR)
   if ( info != 0 ) {
   	  fprintf(stderr," After DGECON INFO = ",info);
   }

   /// call LAPACKE_dgecon with norm I properly (matrix_order=LAPACK_COL_MAJOR)
   if ( info != 0 ) {
   	  fprintf(stderr," After DGECON INFO = ",info);
   }
   c1 = clock();

   printf(" A1norm = %lf, Da1norm = %lf, Rcond = %lf\n", a1norm, da1norm, r1cond); 
   printf(" A2norm = %lf, Da2norm = %lf\n", a2norm, da2norm); 
   printf(" AInorm = %lf, DaInorm = %lf, Rcond = %lf\n", aInorm, daInorm, rIcond); 
   
   printf("Seconds for calculating C.N.: %lf\n",(double)(c1-c0)/(double)CLOCKS_PER_SEC);
  
   return(0);
}