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


// BLAS functions prototypes:

	// Sum of the absolute values.
double dasum_(const int *n, const double *X, const int *incx);

	// DNRM2 returns the euclidean norm of a vector
double dnrm2_(const int *n, const double *X, const int *incx);

/*
   LAPACKE functions prototypes:
   
	! DGETRF computes an LU factorization of a general M-by-N matrix A
    ! using partial pivoting with row interchanges.
    ! The factorization has the form
    !    A = P * L * U
    lapack_int LAPACKE_dgetrf 	( 	int  	matrix_order,
		lapack_int  	m,
		lapack_int  	n,
		double *  	a,
		lapack_int  	lda,
		lapack_int *  	ipiv 
	) 	

	! DGECON estimates the reciprocal of the condition number of a general
    ! real matrix A, in either the 1-norm or the infinity-norm, using
    ! the LU factorization computed by DGETRF.
	lapack_int LAPACKE_dgecon 	( 	int  	matrix_order,
		char  	norm,
		lapack_int  	n,
		const double *  	a,
		lapack_int  	lda,
		double  	anorm,
		double *  	rcond 
	) 	
	
    ! DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or 
	! the largest absolute value of any element of a general rectangular matrix.
	! DLANGE = ( max(abs(A(i,j))), NORM = 'M' or 'm'
    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;
	   // call dasum_ to compute the absolute values of A[][]
	   // use TV for storing the results:
			.	.	.
	   anorm = MAXVAL((*N),TV);
	   free(TV);
	} else if ( Norm[0] == '2' ) {
	   len = (*LDA)*(*N);
	   one = 1;
	   // call dnrm2_ to compute this norm:
			.	.	.
	} else if ( Norm[0] == 'i' || Norm[0] == 'I' ) {
	  // Row wise
	   TV = (double*) malloc(sizeof((double)1.0)*(*LDA));
       if ( TV == NULL ) return(-9);
	   // call dasum_ to compute the absolute values of A[][]
	   // use TV for storing the results:
			.	.	.
	   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;
   
   // use LAPACKE_dgetrf to factorize matrix A (use LAPACK_COL_MAJOR as matrix order):
		.	.	.
   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");  
   
   // use LAPACKE_dgecon to compute rcond with a1norm:
		.	.	.
   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); 

   // use LAPACKE_dgecon to compute rcond with aInorm:
		.	.	.
   if ( info != 0 ) {
   	  fprintf(stderr," After DGECON INFO = ",info);
   }
   printf(" AInorm = %lf, Rcond = %lf\n", aInorm, rcond); 
   
  
   return(0);
}