#include <stdio.h>
#include <stdlib.h>
#include <string.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)

int 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, anorm;
   int j, rc, len, one;
   
   *INFO = 0;
   anorm = 0.0;
   if ( Norm[0] == '1' ) {
	   TV = (double*) malloc(sizeof((double)1.0)*(*N));
       if ( TV == NULL ) return(-9);
	   one = 1;
	   for ( j = 0; j < *N; j++ ) {
	   	   TV[j] = dasum_(LDA, &A[index2D(1,j,(*LDA))],&one);
	   } 
	   anorm = MAXVAL((*N),TV);
	   free(TV);
	} else if ( Norm[0] == '2' ) {
	   len = (*LDA)*(*N);
	   one = 1;
	   anorm = dnrm2_(&len,A,&one);
	} else if ( Norm[0] == 'i' || Norm[0] == 'I' ) {
	  // Row wise
	   TV = (double*) malloc(sizeof((double)1.0)*(*LDA));
       if ( TV == NULL ) return(-9);
	   for ( j = 0; j < (*LDA); j++ ) {
		   TV[j] = dasum_(N, &A[index2D(j,1,(*LDA))], LDA);
	   } 
	   anorm = MAXVAL((*LDA),TV);
	   free(TV);
	} else {
    	*INFO = -1;
    }
	return(anorm);
}


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, rcond, rmin, rmax, v;
   int N, Np1, Nm1, LA, np, pos;
   int i, j, rc;
   lapack_int info;
   
   // 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)
   //
   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);
   } 

//
//  Read A matrix 
//
   for ( i = 0; i < LA; i++ ) {
      for ( j = 0; j < N; j++ ) {
   	     fscanf(stdin,"%lf ",&A[index2D(j,i,LA)]);
   	  }
   }

   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..."); 
   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);
   }
   
   printf("Computing matrix factorization...\n");   
   for ( i = 0; i < LA; i++ ) Ipiv[i] = 0.0;
   
   info = LAPACKE_dgetrf(LAPACK_COL_MAJOR, LA, LA, A, LA, Ipiv );
   if ( info != 0 ) {
   	  fprintf(stderr," After DGETRF INFO = %d\n",info);
   }
   
   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");  
   
   info = LAPACKE_dgecon(LAPACK_COL_MAJOR,'1',LA,A,LA,a1norm, &rcond);	
   if ( info != 0 ) {
   	  fprintf(stderr," After DGECON INFO = ",info);
   }
   printf(" A1norm = %lf, Rcond = %lf\n", a1norm, rcond); 
   
   if ( info != 0 ) {
   	  fprintf(stderr," After DGECON INFO = ",info);
   }
   printf(" A2norm = %lf\n", a2norm); 

   info = LAPACKE_dgecon(LAPACK_COL_MAJOR,'I',LA,A,LA,aInorm, &rcond);	
   if ( info != 0 ) {
   	  fprintf(stderr," After DGECON INFO = ",info);
   }
   printf(" AInorm = %lf, Rcond = %lf\n", aInorm, rcond); 
   
  
   return(0);
}