JCMTReduction/linear.cpp at master · hythloda/JCMTReduction · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#include <string.h>
#include "linear.h"
#include <stdio.h>
#include <stdarg.h>

Matrix Matrix::inv()
{
  static Matrix Idt;
  static Matrix Mtrx;
  unsigned short s,r,c;
  double fraction;

  //create inverse matrix
  for (r=0; r<N; r++)
    for (c=0; c<N; c++)
    {
      Idt.pt[r][c] = r == c;
      Mtrx[r][c] = pt[r][c];
    }

// run the matrix inversion algorithm
  for (s=0; s<N; s++)
    for (r=0; r<N; r++)
      if (r!=s)
	for (fraction = Mtrx[r][s]/Mtrx[s][s], c=0; c<N; c++)
	  {
            Mtrx[r][c] =    Mtrx[r][c] -  Mtrx[s][c]*fraction;
	    Idt[r][c] =   Idt[r][c] - Idt[s][c]*fraction;
	  }


  for (r=0; r<N; r++)
    for (c=0; c<N; c++)
     {
       //           printf("Idt[%d][%d] = %lf\n", r,c,Mtrx[r][c]/Mtrx[r][r]);
             Idt[r][c] = Idt[r][c]/Mtrx[r][r];
     }

  return Idt;

}


char *Vector::print()
{
  static char string[200];
  char buffer[30];

  for (unsigned short i=0; i<N; i++)
    {
      sprintf(buffer,"%f",pt[i]);
      strcat(string,buffer);
      strcat(string,"  ");
    }

  return string;
}


Vector vector(double zeroeth,...)
{
  Vector result;
  va_list list;

  result[0] = zeroeth;
  va_start(list, zeroeth);
  for(unsigned short i=1; i<N; i++)
      result[i] = va_arg(list, double);
  va_end(list);

  return result;
}