#include<math.h>
#include<stdio.h>
#include<time.h>
#include<stdlib.h>
main(){
       double X[200][200],Q[200][200],alfa[200][200],z[200],norm,prvok,rez,orto;
       double rez1, orto1, Q1[200][200], ALFA1[200][200];
       int n,i,j,k;
       //vstup = rad matice
       printf("Zadaj rozmer matice: ");
       scanf("%d",&n);
       //nahodna matica
       //start generatora nahodnych cisel
       srand((unsigned)time(NULL));
       for(j=0;j<n;j++){
          //cyklus j=stlpec
          for(i=0;i<n;i++){
             X[i][j]=(-0.5+rand()/(double)RAND_MAX)*10.E+1;
             Q[i][j]=(double)0.;
             alfa[i][j]=(double)0.;
             Q1[i][j]=(double)0.;
             ALFA1[i][j]=(double)0.;
          }//i
       }//j
//**********************************
       for(j=0;j<n;j++){
          //cyklus j = stlpec
                        for(i=0;i<j;i++){
                                         //skalarny sucin qi*xj
                                         alfa[i][j]=(double)0.;
                                         for(k=0;k<n;k++)
                                         alfa[i][j]+=Q[k][i]*X[k][j];
                                         }
                                         //i
                                         //pomocny vektor z
                                         for(k=0;k<n;k++)
                                         {
                                                          z[k]=(double)0.;
                                                          for(i=0;i<j;i++)
                                                          z[k]+=alfa[i][j]*Q[k][i];
                                                          }
                                                          //k
                                                          //ortogonalizacia
                                                          for(k=0;k<n;k++)Q[k][j]=X[k][j]-z[k];
                                                          //norma Q[*][j]
                                                          norm=(double)0.;
                                                          for(k=0;k<n;k++)norm+=Q[k][j]*Q[k][j];
                                                          norm=sqrt(norm);
                                                          alfa[j][j]=norm;
                                                          //normovanie vektora Q[*][j]
                                                          for(k=0;k<n;k++)Q[k][j]=Q[k][j]/norm;
                                                          }//j
                                                          //strata ortogonality
                                                          orto=(double)0.;
                                                          for(i=0;i<n;i++){
                                                                           for(j=0;j<n;j++){
                                                                                            prvok=(double)0.;
                                                                                            for(k=0;k<n;k++)prvok+=Q[k][i]*Q[k][j];
                                                                                            if(i==j)prvok=1.-prvok;
                                                                                            else prvok= -prvok;
                                                                                            orto+=prvok*prvok;
                                                                                            }//j
                                                                                            }//i
                                                                                            orto=sqrt(orto);
                                                                                            
  //presnost QR fakt.=Frob.norma rezidua (X - Q*ALFA)
  rez=(double)0.;
  for(i=0; i<n ;i++){
      for(j=0; j<n; j++){
         prvok = (double)0.;
         for (k=0; k<=j; k++) prvok += Q[i][k] * alfa[k][j];
         prvok = X[i][j] - prvok;
         rez += prvok * prvok;
      }
  }
  rez = sqrt(rez);
  
  //Modifikovany Gramm-Schmidt
  for(j=0; j<n; j++){
     //cyklus cez stlpce

     //kopia X[*][j[ do Q1[*][j]
     for(k=0; k<n; k++) Q1[k][j] = X[k][j];
     
     //skalarny sucin
     for (i=0; i<j; i++){
         ALFA1[i][j] = (double)0.;
         for(k=0; k<n; k++) ALFA1[i][j] += Q1[k][i] * Q1[k][j];
         
         //hned pouzijem ALFA1[i][j] na ortogonalizaciu
         for(k=0; k<n; k++) Q1[k][j] = Q1[k][j] - ALFA1[i][j] * Q1[k][i];
     }//koniec cyklu cez 1
     
     //norma Q1[*][j]
     norm = (double)0.;
     for(k=0; k<n; k++) norm += Q1[k][j]*Q1[k][j];
     
     norm=sqrt(norm);
     
     ALFA1[j][j] = norm; //diag. prvok hornej torj. ALFA1
     
     //normovanie vektora
     for(k=0; k<n; k++) Q1[k][j] = Q1[k][j] / norm;
  }//cyklus j
  
  //strata ortogonality = Frobeniova norma matice (I-Q1^T*Q1)

  orto1=(double)0.;
  for(i=0;i<n;i++){
     for(j=0;j<n;j++){
        prvok=(double)0.;
        for(k=0;k<n;k++)prvok+=Q1[k][i]*Q1[k][j];
        if(i==j)prvok=1.-prvok;
        else prvok= -prvok;
        orto1+=prvok*prvok;
     }//j
  }//i
  orto1=sqrt(orto1);

  //presnost QR fakt.=Frob.norma rezidua (X - Q1*ALFA1)
  rez1=(double)0.;
  for(i=0; i<n ;i++){
      for(j=0; j<n; j++){
         prvok = (double)0.;
         for (k=0; k<=j; k++) prvok += Q1[i][k] * ALFA1[k][j];
         prvok = X[i][j] - prvok;
         rez1 += prvok * prvok;
      }
  }
  rez1=sqrt(rez1);
  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  
  
  //vypis straty ortogonality a rezidua
  printf("\n\nCGS: strata OG     : %G \n", orto);
  printf("CGS: reziduum     : %G \n", rez);
  
  printf("MGS: strata OG    : %G \n", orto1);
  printf("MGS: reziduum     : %G \n", rez1);

  system("pause");
}