/* FILE: matutl6.c begins */

/* Muutettu 20.1.2000

/* rdm, ranmat have been converted from double to float */
/* Matriisin k\"asittelyfunktioita kurssille
   Numeeriset menetelm\"at ja C-kieli NRC00
  (matriisin tyyppi float **, kuten tiedostossa nrutil.c) 
   Kutsuvassa ohjelmassa pit\"a\"a olla m\"a\"arittelyt:  
#include "nr.h"
#include "nrutil.h"
#include "nrutil.c"
#include "pythag.c"
#include "svdcmp.c"
*/

#ifndef _NR_H_
#define _NR_H_
#include "../util/nr.h"
#endif

#ifndef _NR_UTILS_H_
#define _NR_UTILS_H_
#include "../util/nrutil.h"
#endif



#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include <time.h>
#include <math.h>


char *getfname(char *fname)
{
  FILE *fp;
  int count=1;
  char *name,*beg,*end,*point;
  
  name=(char *) malloc(20);
  beg=(char *) malloc(20);
  sprintf(name,"%s",fname);
  sprintf(beg,"%s",fname);
  point=strchr(beg,'.');
  if(point!=NULL)
    *point='\0';
  end=".dat";
  while(!((fp=fopen(name,"r"))==NULL)){ 
    sprintf(name,"%s%d%s",beg,count++,end);
    fclose(fp);
  }
  return name;
}


/* This routine makes the seed for the random generator */
void init_srand(void)
{
  time_t *tp, aika;
  unsigned int siemen;	
  
  tp =  (time_t *) malloc(sizeof(time_t));
  aika = time(tp); 
  siemen =  (unsigned int) aika%10000;
  srand(siemen);
  return;
}

/*
double rdm(double y1, double y2)
{
  double x,der;
  x=rand();
  der=(double) (y2 - y1)/RAND_MAX;
  return (der*x + y1);
}
*/

float rdm(float y1, float y2)
{
  float x,der;
  x=rand();
  der=(float) (y2 - y1)/RAND_MAX;
  return (der*x + y1);
}


/* Returns random matrix of given size */
float **ranmat(long m,long n)
{
  int i,j;	
  float **c;
  float rd1=-10.,rd2=10.; /* Was: double */
  c=matrix(1,m,1,n); 
  for (i=1;i<=m;i++){
    for (j=1;j<=n;j++){
      c[i][j]= rdm(rd1,rd2);}}
  return c;
}


/* Returns random matrix of given size */
/* float **ranmat2(long m,long n, double rd1, double rd2) */
float **ranmat2(long m,long n, float rd1, float rd2)
{
  int i,j;	
  float **c;
  c=matrix(1,m,1,n); 
  for (i=1;i<=m;i++){
    for (j=1;j<=n;j++){
      c[i][j]= rdm(rd1,rd2);}}
  return c;
}


/* Displays a matrix on the screen. */
void showmat(float **a,long m,long n)
{
  int i,j;
  printf("\n");
  for (i=1;i<=m;i++){
    for (j=1;j<=n;j++)
      printf("%- 12.6f",a[i][j]);
    printf("\n");
  }
}

/* Prompts the user to enter a matrix entry by entry */
float **entermat(long *mm,long *nn)
{
  int i,j;
  long m,n;
  float **a;
  FILE *fp;
  printf("Enter integers m and n (<8) ");
  scanf(" %ld %ld",&m,&n);
  a=matrix(1,m,1,n);
  for (i=1;i<=m;i++){
    for (j=1;j<=n;j++){
      printf("a[%d,%d] : ",i,j);
      scanf(" %f",&a[i][j]);}}
  *mm=m;
  *nn=n;
  return a;
}

/* Writes a matrix into a file */
char *putmat(float **a,long m,long n,char *tied)
{
  FILE *fp;
  int i,j;
  char *fname;
  
  fname=getfname(tied);
  if ((fp=fopen(fname,"w"))==NULL){
    fprintf(stderr,"Cannot open file %s\n",fname);
    exit(1);
  }
  fprintf(fp,"%ld %ld\n",m,n);
  for (i=1;i<=m;i++){
    for (j=1;j<=n;j++) fprintf(fp,"%f ",a[i][j]);
    fprintf(fp,"\n");}
  fclose(fp);
  return fname;
}

/* Reads a matrix from a file. */
float **getmat(long *mm,long *nn,char *tied)
{
  FILE *fp;
  float **a;
  int i,j;
  long m,n;
  char dum[80];
  if ((fp=fopen(tied,"r"))==NULL){
    fprintf(stderr,"Cannot open file %s\n",tied);
    exit(1);
  } 
  if (fscanf(fp,"%ld %ld",&m,&n)!=2) 
    printf("Bad matrix file");
  a=matrix(1,m,1,n);
  for (i=1;i<=m;i++){
    for (j=1;j<=n;j++) 
      fscanf(fp," %f",&a[i][j]);
    fscanf(fp,"\n");}
  fclose(fp);
  *mm=m;*nn=n;
  return a;
}

/* Multiplies two matrices */
float **matmul(float **a,float **b,long n,long p,long m)
{
  float **c;
  int i,j,k;
  c=matrix(1,n,1,m);
  for (i=1;i<=n;i++){
    for (j=1;j<=m;j++){
      c[i][j]=0;
      for (k=1;k<=p;k++) 
	c[i][j]+=(a[i][k]*b[k][j]);
    }
  }
  return c;
}

/* Matrix Norm: Sum |a(i,j)| */
float mnorm1(float **a,long m, long n)
{
  float s=0.0;
  int i,j;
  for (i=1;i<=m;i++){
    for (j=1;j<=n;j++){
      s=s+fabs(a[i][j]);}}
  return s;
}



/* Returns matrix with all elements = 1, of given size */
float **ones(long m,long n)
{
  int i,j;	
  float **c;
  c=matrix(1,m,1,n); 
  for (i=1;i<=m;i++){
    for (j=1;j<=n;j++){
      c[i][j]= 1.0;}}
  return c;
}


/* Returns matrix with element =1 on diagonal,
   0 elsewhere, of given size */
float **unitmat(long m,long n)
{
  int i,j;	
  float **c;
  c=matrix(1,m,1,n); 
  for (i=1;i<=m;i++){
    for (j=1;j<=n;j++){
      if (i==j) {c[i][j] =1.0;}	
      else c[i][j]= 0.0;}}
  return c;
}

/* Transpoose of a matrix */
float **transp(float **a,long m,long n)
{
  float **c;
  int i,j,k;
  c=matrix(1,n,1,m);
  for (i=1;i<=m;i++) 
    for (j=1;j<=n;j++) c[j][i]=a[i][j];
  return c;
}

double cond(float **aa, long m, long n)
     /* Returns condition number of n by n matrix */
     /* In the calling program must define:
	#include "pythag.c"
	#include "svdcmp.c"
     */
{
  int i,j;
  float *w,**u,**v;
  double f,wmin, wmax ;
  w=vector(1,n);
  u=matrix(1,m,1,n);
  v=matrix(1,n,1,n);
  if (m !=n) {
    printf("cond not applicable, m= %d, n= %d \n",m,n);
    exit(1);}
  for (i=1;i<=m;i++)
    for (j=1;j<=n;j++) {                    
      u[i][j]=aa[i][j];}
  svdcmp(u,m,n,w,v);
  wmax=w[1];  wmin= w[1];
  for (i=1;i<=n;i++){
    if (w[i] >wmax) {wmax =w[i];}        
    if (w[i] <wmin){ wmin =w[i];}
  }
  if (wmin <1e-20) wmin =1e-20;
  free_matrix(v,1,n,1,n);
  free_matrix(u,1,m,1,n);
  free_vector(w,1,n);
  return(wmax/wmin);
} 
/* FILE: matutl6.c ends */