matrixjpl.c

计量工具箱

/*===================================================================
 matrixjpl.c

 carried out linear algebra routines

===================================================================*/

#include <math.h>
#include <malloc.h>
#include <stddef.h>
#include <stdio.h>
#include "matrixjpl.h"

double dabs(double in)
{
	double out;

	if (in > 0.0)
		out = in;
    if (in < 0.0)
		 out = -in;
    if (in == 0.0)
        out = 0.0;
  
return out;
}


void vextract(int n, int col, double **mat, double *vec)
{ // extracts column 'col' from a double matrix
	int i;
	// N.B. col is 0-indexed as in C

	for(i=0; i<n; i++)
		vec[i] = mat[i][col];
}

void vinsert(int n, int col, double *vec, double **mat)
{ // inserts vec into 'col' of a double matrix
	int i;
	// N.B. col is 0-indexed as in C

	for(i=0; i<n; i++)
		mat[i][col] = vec[i];
}

void ivextract(int n, int col, int **mat, int *vec)
{ // extracts column 'col' from an integer matrix
	int i;
	// N.B. col is 0-indexed as in C

	for(i=0; i<n; i++)
		vec[i] = mat[i][col];
}


void vcopy(int n, double *vec1, double *vec2)
{ // copies vec1 to vec2
	int i;

	for(i=0; i<n; i++)
		vec2[i] = vec1[i];
}




// vprint prints double vectors
void vprint(int n, double *vec)
{
int i;

for(i=0; i<n; i++)
printf(" %8.3lf \n",vec[i]);

printf("\n");

}

// ivprint prints integer vectors
void ivprint(int n, int *vec)
{
int i;

for(i=0; i<n; i++)
printf(" %8d \n",vec[i]);

printf("\n");

}


// mprint prints double matrices
void mprint(int n, int m, double **vec)
{
int i, j;

for(i=0; i<n; i++){
	for(j=0; j<m; j++)
	printf(" %8.3lf ",vec[i][j]);
	printf(" \n");
}

printf("\n");


}

// imprint prints integer matrices
void imprint(int n, int m, int **vec)
{
int i, j;

for(i=0; i<n; i++){
	for(j=0; j<m; j++)
	printf(" %8d ",vec[i][j]);
	printf(" \n");
}

printf("\n");

}

// permute() and permute_inverse() perform
/* In-place Permutations 

   permute:    OUT[i]       = IN[perm[i]]     i = 0 .. N-1
   invpermute: OUT[perm[i]] = IN[i]           i = 0 .. N-1

   PERM is an index map, i.e. a vector which contains a permutation of
   the integers 0 .. N-1.

   From Knuth "Sorting and Searching", Volume 3 (3rd ed), Section 5.2
   Exercise 10 (answers), p 617
*/

void permute(int n, int *p, double *data)
{
  int i, k, pk;
  double r1, t;
  
  for (i = 0; i < n; i++)
    {
      k = p[i];
      
      while (k > i) 
        k = p[k];
      
      if (k < i)
        continue ;
      
      /* Now have k == i, i.e the least in its cycle */
      
      pk = p[k];
      
      if (pk == i)
        continue ;
      
      /* shuffle the elements of the cycle */
      
      t = data[i];
      
      while (pk != i)
        {
          r1 = data[pk];
          data[k] = r1;
          k = pk;
          pk = p[k];
        };
      
      data[k] = t;
    }
}

void permute_inverse (int n, int *p, double *data)
{
  int i, k, pk;

  double r1, t;

  for (i = 0; i < n; i++)
    {
      k = p[i];
          
      while (k > i) 
        k = p[k];

      if (k < i)
        continue ;
      
      /* Now have k == i, i.e the least in its cycle */

      pk = p[k];

      if (pk == i)
        continue ;
      
      /* shuffle the elements of the cycle in the inverse direction */
      
      t = data[k];

      while (pk != i)
        {
          r1 = data[pk];
          data[pk] = t;
          k = pk;
          pk = p[k];
          t = r1;
        };
      
      data[pk] = t;
    }

}




void matmulc(double z[], double x[], double y[], int size,
             int m1, int m2, int n1, int n2)
{// mutliplies non-conformable matrices so long as they
 // are conformable in a least one dimension
  int i, j;
  int flag = 0;
  
  /* case of two equal matrices */
  if (m1 == m2  &&  n1 == n2){
  for (i=0; i <= size-1; i++){
  z[i] = x[i]*y[i];
  }
    flag = 1;
  }
  if (m1 == m2 && n2 == 1) {
  size = 0;
  for (j=0; j<= n1-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = x[size] * y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (m1 == m2 && n1 == 1) {
  size = 0;
  for (j=0; j<= n2-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = y[size] * x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (n1 == n2 && m2 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m1-1; j++){
  z[size] = x[size] * y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
  if (n1 == n2 && m1 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m2-1; j++){
  z[size] = y[size] * x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
   if (flag == 0){
      printf("matmulc: nonconformability of matrices \n");
   }

}

void matdivc(double z[], double x[], double y[], int size,
             int m1, int m2, int n1, int n2)
{// divides non-conformable matrices so long as they
 // are conformable in a least one dimension
  int i, j;
  int flag = 0;
  
  /* case of two equal matrices */
  if (m1 == m2  &&  n1 == n2){
  for (i=0; i <= size-1; i++){
  z[i] = x[i]/y[i];
  }
    flag = 1;
  }
  if (m1 == m2 && n2 == 1) {
  size = 0;
  for (j=0; j<= n1-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = x[size] / y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (m1 == m2 && n1 == 1) {
  size = 0;
  for (j=0; j<= n2-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = y[size] / x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (n1 == n2 && m2 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m1-1; j++){
  z[size] = x[size] / y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
  if (n1 == n2 && m1 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m2-1; j++){
  z[size] = y[size] / x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
   if (flag == 0){
      printf("matdivc: nonconformability of matrices \n");
   }

}

void mataddc(double z[], double x[], double y[], int size,
             int m1, int m2, int n1, int n2)
{// adds non-conformable matrices so long as they
 // are conformable in a least one dimension
  int i, j;
  int flag = 0;
  
  /* case of two equal matrices */
  if (m1 == m2  &&  n1 == n2){
  for (i=0; i <= size-1; i++){
  z[i] = x[i]+y[i];
  }
    flag = 1;
  }
  if (m1 == m2 && n2 == 1) {
  size = 0;
  for (j=0; j<= n1-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = x[size] + y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (m1 == m2 && n1 == 1) {
  size = 0;
  for (j=0; j<= n2-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = y[size] + x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (n1 == n2 && m2 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m1-1; j++){
  z[size] = x[size] + y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
  if (n1 == n2 && m1 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m2-1; j++){
  z[size] = y[size] + x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
   if (flag == 0){
      printf("mataddc: nonconformability of matrices \n");
   }

}

void matsubc(double z[], double x[], double y[], int size,
             int m1, int m2, int n1, int n2)
{// subtracts non-conformable matrices so long as they
 // are conformable in a least one dimension
  int i, j;
  int flag = 0;
  
  /* case of two equal matrices */
  if (m1 == m2  &&  n1 == n2){
  for (i=0; i <= size-1; i++){
  z[i] = x[i] - y[i];
  }
    flag = 1;
  }
  if (m1 == m2 && n2 == 1) {
  size = 0;
  for (j=0; j<= n1-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = x[size] - y[i];
  size = size + 1;