/* Najprej definiramo stevilo korakov, v katerih naj racuna
   po formuli. 2.000 je "standardno" stevilo korakov za test hitrosti.

   DIMEN doloca dimenzijo matrike; to se ne sme spreminjat, ker je zacetna
   matrika (katere inverz racunamo) staticno definirana!
*/

#define N_ITER 2000
#define DIMEN 16

/* Potem nekaj pogojev ...
   Ce je definiran DEBUG, potem sproti izpisuje kaj pocne.

   Ce je definiran WANT_RESULTS, potem na koncu izpise rezultate.

   VAX doloca, da prevajamo na VAX/VMS masini in ima drugace dolocen
   korak sistemskega timerja (druga simbolicna spremenljivka)

   RAW doloca, kako naj se dostopa do matrike; ali je lahko v enem
   ciklu dostop le do ene vrstice (v tem primeru rabimo zacasno
   spremenljivko, v katero spravimo eno vrstico) ali pa lahko v eni
   zanki dostopamo do dveh vrstic hkrati.
*/

/* #define DEBUG */
/* #define WANT_RESULTS */
/* #define VAX */
#define RAW


/* No, sedaj pa zares.... */

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

main()
{
   /* v matriki cc je pravzaprav matrika, katere inverz racunamo */
   double cc[DIMEN][DIMEN]={
       { 1.0, 0.2,-0.5,1.25, 0.0,-0.3, 1.5,-1.5,1.25,-2.0, 0.0,-0.7, 1.1,-0.1, 0.0, 1.1},
       { 2.3,-1.6, 0.7,-0.8, 0.2,0.75,-0.3, 0.0, 1.0,-0.9, 1.1, 0.1,-0.3,-0.2, 0.3,-0.6},
       { 0.0, 1.6,-1.3, 0.7, 0.0, 1.2,-0.5, 0.6, 1.2,-0.5, 0.4,-0.3, 0.7, 0.0, 0.0, 0.0},
       { 1.1, 1.0,-0.3, 1.1, 0.2, 0.0, 0.0,-0.2, 0.0,-0.3,-0.6, 0.0, 1.3, 0.0, 1.1, 0.0},
       {-0.3, 0.1,-0.2,-2.1, 1.8, 0.0,-0.3,-0.1, 0.0, 0.0, 1.1, 0.0, 1.3,-0.9, 0.2, 0.7},
       { 1.0, 0.0,-0.4, 0.0, 0.6,-0.3, 0.0, 0.0, 1.0,-0.2, 0.0, 0.0, 1.5,-1.4, 0.0, 0.3},
       { 0.0,-0.2, 1.4, 0.2,-0.3, 0.0, 1.7, 0.1,-0.8, 1.1, 0.0, 0.1,-1.3,-1.2, 0.2,-0.7},
       { 0.0, 0.0, 1.3,-1.9,-0.2, 0.7, 0.0, 1.2,-0.6, 0.5, 0.1, 0.0, 1.2,-1.3, 0.0, 0.0},
       { 1.0, 0.0, 0.0, 0.0, 0.0, 0.1, 0.3,-0.4, 0.2, 1.3,-1.5, 0.2, 0.0,-0.1, 1.2, 1.2},
       { 0.0, 0.1, 0.2, 0.3, 0.0,-0.4,-0.5, 0.0, 0.0, 0.0, 0.0, 1.9, 0.3,-0.3, 0.0, 0.0},
       {-0.2, 1.7, 0.7, 0.0,-0.1, 0.4, 0.0, 0.0, 1.2,-0.7, 0.1, 0.0, 0.0,-0.1, 1.1,-0.4},
       { 0.0, 1.1, 1.0,-0.9, 0.8,-0.7, 0.2,-0.1, 0.0, 0.3, 0.2, 0.7,-0.3, 0.0, 1.1,-1.9},
       { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.1, 0.0, 1.2, 0.1, 0.2,-0.3,-0.4,-0.5, 0.6, 0.7},
       { 0.1, 0.2, 0.3, 1.2,-0.6,-0.5,-1.3, 0.1,-0.6,-0.6,-0.7, 0.0,-1.2, 2.0, 0.0,-1.7},
       { 0.0, 0.0, 1.1, 0.2,-0.3,-0.3, 0.4, 0.5, 0.2, 1.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1},
       { 0.1, 0.1, 0.0, 0.0,-0.2, 1.7, 1.6,-0.5,-1.2,-0.3, 0.1, 0.0, 0.0, 1.3, 0.2,-1.2}
   };
   double ab[DIMEN][2*DIMEN];   /* matriki a in b (inverz a) zdruzeni v eno */

#ifdef RAW
   double vrsta[2*DIMEN];       /* pomozna spremenljivka za eno vrsto, */
                                /* rabimo jo samo v primeru RAW */
#endif

   double diag,tmp;       /* pomozna speremenljivka */

   int i,j,k;
   long n;
   clock_t start,end;



#ifdef DEBUG
   printf ("Iteracija:  vrstica:\n");
#endif

   start=clock();
   for (n=0;n<N_ITER;n++)
   {
#ifdef DEBUG
      printf ("   %5d\r",n);
#endif
      for (i=0;i<DIMEN;i++)
         for (j=0;j<DIMEN;j++)
            ab[i][j]=cc[i][j];

      for (i=0;i<DIMEN;i++)
      {
         for (j=0;j<DIMEN;j++)
            ab[i][DIMEN+j]=0.0;
         ab[i][DIMEN+i]=1.0;
      }
      for (j=0;j<DIMEN;j++)
      {
#ifdef DEBUG
         printf ("   %5d  %3d\r",n,j);
#endif
         diag=ab[j][j];
         for (i=0;i<2*DIMEN;i++)
            ab[j][i]=ab[j][i]/diag;
         for (i=0;i<j;i++)           /* obdelavo matrike razstavim na vrstice */
                                     /* pred trenutno in tiste za, trenutno */
                                     /* moramo pustiti pri miru */
         {
            tmp=ab[i][j];
#ifdef RAW
            for (k=0;k<2*DIMEN;k++)
               vrsta[k]=ab[j][k]*tmp;

            for (k=0;k<2*DIMEN;k++)
               ab[i][k]=ab[i][k]-vrsta[k];
#else
            for (k=0;k<2*DIMEN;k++)
               ab[i][k]=ab[i][k]-ab[j][k]*tmp;
#endif
         }

         for (i=j+1;i<DIMEN;i++)
         {
            tmp=ab[i][j];
#ifdef RAW
            for (k=0;k<2*DIMEN;k++)
               vrsta[k]=ab[j][k]*tmp;

            for (k=0;k<2*DIMEN;k++)
               ab[i][k]=ab[i][k]-vrsta[k];
#else
            for (k=0;k<2*DIMEN;k++)
               ab[i][k]=ab[i][k]-ab[j][k]*tmp;
#endif
         }
      }
   }
   end=clock();

#ifdef DEBUG
   printf ("\n");
#endif

#ifdef WANT_RESULTS
   printf ("Matrika:\n\n");
   for (j=0;j<DIMEN;j++)
   {
      for (i=0;i<DIMEN;i++)
         printf ("%5.3f  ",cc[j][i]);
      printf ("\n");
   }
   printf ("\n");

   printf ("Inverz matrike:\n\n");
   for (j=0;j<DIMEN;j++)
   {
      for (i=DIMEN;i<2*DIMEN;i++)
         printf ("%5.3f  ",ab[j][i]);
      printf ("\n");
   }
   printf ("\n");
#endif

#ifdef VAX
   printf ("  Porabljen cas: %8.4fs\n",((double)end-start)/CLK_TCK);
#else
   printf ("  Porabljen cas: %8.4fs\n",((double)end-start)/CLOCKS_PER_SEC);
#endif

   return (0);
      
}