prova.c

Go to the documentation of this file.

/*
    This file is part of the Ansel project.
    Copyright (C) 2023 Davide Patria.
 
    Ansel is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    Ansel is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
    along with Ansel.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#include "polar_decomposition.h"
#include <stdio.h>
 
void print_float(TYPE num);
 
void check_orthonormal(TYPE Q[3][3]);
 
void stampa_3x3(TYPE Q[3][3]);
 
void copy(TYPE *dest, TYPE *source, size_t num_el);
 
#ifndef TEST_N
#define TEST_N 1
#endif
 
int main()
{
 
  TYPE A[3][3];
  // simple case one
  switch(TEST_N)
  {
    case 1:
    {
      TYPE temp[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
      copy(*A, *temp, 9);
    }
    case 2:
    {
      TYPE temp[3][3] = { { 0.1, 0.1, 0.3 }, { 0.2, 0.1, 0.2 }, { 0.3, 0.0, 0.1 } };
      copy(*A, *temp, 9);
    }
    // TODO: check and use here the maledict test mentioned in the cpp version
    case 3:
    {
      TYPE temp[3][3] = { { -1, 2, 0 }, { 1.3, 0.4, 3.4 }, { 0.3, 2.4, 3.4 } };
      copy(*A, *temp, 9);
    }
  }
 
  TYPE Q[3][3];
  TYPE H[3][3];
 
  // stampa_3x3(A);
  polar_decomposition(A, Q, H);
 
  // stampa_3x3(A);
  stampa_3x3(Q);
  stampa_3x3(H);
 
  matrix_multiply(*(Q), *(H), *(A), 3, 3, 3);
 
  // stampa_3x3(A);
 
  check_orthonormal(Q);
 
  return 0;
}
 
void stampa_3x3(TYPE Q[3][3])
{
  for(size_t riga = 0; riga < 3; riga++)
  {
    for(size_t col = 0; col < 3; col++)
    {
      printf("%f ", Q[riga][col]);
      if(col == 2) printf("\n");
    }
  }
  printf("\n");
}
 
void check_orthonormal(TYPE Q[3][3])
{
  TYPE scalar01 = 0;
  TYPE scalar12 = 0;
  TYPE scalar02 = 0;
 
  for(size_t i = 0; i < 3; i++)
  {
    scalar01 += Q[0][i] * Q[1][i];
    scalar12 += Q[1][i] * Q[2][i];
    scalar02 += Q[0][i] * Q[2][i];
  }
  print_float(scalar01);
  print_float(scalar12);
  print_float(scalar02);
}
 
void print_float(TYPE num)
{
  printf("scalar12: %f\n", num);
}