/*
    This file is part of the FElt finite element analysis package.
    Copyright (C) 1993-2000 Jason I. Gobat and Darren C. Atkinson

    This program 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 2 of the License, or
    (at your option) any later version.

    This program 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 this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/

/************************************************************************
 * File:	data.c	
 *	
 * Description:	
 *		
 ************************************************************************/

# include <stdio.h>
# include <stdlib.h>
# include "matrix.h"
# include "error.h"

double mdata (A, row, col)
   Matrix	A;
   unsigned	row;
   unsigned	col;
{
   unsigned	height;
   unsigned	temp;

   if (IsFull (A))
      return A -> data [row][col];
   else if (IsCompact (A)) {
      if (row > col) {
         temp = col;
         col = row;
         row = temp;
      }

      if (col == 1)
         height = 1;
      else
         height = A -> diag [col] - A -> diag [col - 1];

      if (row > col - height)
         return A -> data [A -> diag [col] + row - col][1]; 
      else
         return 0;
   }

   return 0;
}

Matrix CreateSubsectionMatrix (a, sr, sc, er, ec)
   Matrix	a;
   unsigned	sr, sc;
   unsigned	er, ec;
{
   Matrix	b;
   unsigned	i;
   Matrix	root;

   if (er > a -> nrows || ec > a -> ncols || sc < 1 || sr < 1)
      return NULL;

   if (IsCompact (a))
      return NULL;

   b = (struct matrix *) malloc (sizeof (struct matrix));
   if (b == NULL)
	Fatal ("unable to allocate subsection matrix");

   b -> nrows = er - sr + 1;
   b -> ncols = ec - sc + 1;

   b -> data = (double **) malloc (sizeof (double *) * b -> nrows);
   if (b -> data == NULL)
	Fatal ("unable to allocate subsection matrix");
   b -> data --;

   for (i = 1 ; i <= b -> nrows ; i++)
      b -> data [i] = &(a -> data [i + sr - 1][sc - 1]); 

   b -> refcount = 0;
   b -> size = 0;

   root = a;
   while (root -> parent)
      root = root -> parent;

   b -> parent = root;
   root -> refcount ++;

   return b;
}

Matrix CreateFullMatrix (rows, cols)
   unsigned	rows;
   unsigned	cols;
{
   unsigned	i;
   Matrix	m;

   m = (struct matrix *) malloc (sizeof (struct matrix));
   if (m == NULL) {
        fprintf (stderr,"failure point 1: r = %d, c = %d\n", rows, cols);
	Fatal ("unable to allocate full matrix");
   }

   m -> data = (double **) malloc (sizeof (double *) * rows);
   if (m -> data == NULL) {
        fprintf (stderr,"failure point 2: r = %d, c = %d\n", rows, cols);
	Fatal ("unable to allocate full matrix");
   }

   m -> data --;

   m -> data [1] = (double *) malloc (sizeof (double) * rows * cols);
   if (m -> data [1] == NULL) {
        fprintf (stderr,"failure point 3: r = %d, c = %d\n", rows, cols);
	Fatal ("unable to allocate full matrix");
   }

   m -> data [1] --;
   for (i = 2 ; i <= rows ; i++)
      m -> data [i] = m -> data [i-1] + cols;

   m -> nrows = rows;
   m -> ncols = cols;
   m -> size = 0;
   m -> refcount = 1;
   m -> parent = NULL;

   return m;
}

Matrix CreateRowVector (size)
   unsigned	size;
{
   return CreateFullMatrix (1, size);
}

Matrix CreateColumnVector (size)
   unsigned	size;
{
   return CreateFullMatrix (size, 1);
}

void DestroyMatrix (m)
   Matrix	m;
{
   if (m -> parent != NULL) {
      m -> parent -> refcount --;
 
      if (m -> parent -> refcount == 0)
         DestroyMatrix (m -> parent);
      
      m -> data ++;
      free (m -> data);
      free (m); 

      return;
   } 
   else if (-- m -> refcount)
      return;

   m -> data [1] ++;
   free (m -> data [1]);

   m -> data ++;
   free (m -> data);

   if (IsCompact (m)) {
      m -> diag ++;
      free (m -> diag);
   }
   
   free (m);
}

Matrix CreateCompactMatrix (rows, cols, size, diag)
   unsigned	rows;
   unsigned 	cols;
   unsigned	size;
   unsigned	*diag;
{
   Matrix	A;

   A = CreateFullMatrix (size, 1);

   A -> nrows = rows; 
   A -> ncols = cols;
   A -> size = size;
   A -> parent = NULL;

   if (diag != NULL) {
      A -> diag = (unsigned *) malloc (sizeof (unsigned) * rows);
      if (A -> diag == NULL)
         Fatal ("unable to create compact matrix");

      A -> diag --;

      memcpy ((void *) &(A -> diag[1]), (void *) &(diag[1]), 
              sizeof(unsigned)*rows);
   }

   return A;
}

Matrix CreateCopyMatrix (a)
   Matrix	a;
{
   Matrix	b;
   unsigned	size;
   unsigned	rows, cols;

   rows = a -> nrows;
   cols = a -> ncols;

   if (IsCompact(a)) {
      size = Msize(a)*sizeof(double);
      b = CreateCompactMatrix (rows, cols, Msize(a), NULL); 
      b -> diag = (unsigned *) malloc (sizeof(unsigned) * rows);
      if (b -> diag == NULL)
	Fatal ("unable to create compact matrix");
      b -> diag --;
      memcpy ((void *) &(b -> diag[1]), (void *) &(a -> diag[1]), 
              sizeof(unsigned)*rows);
   }
   else { 
      size = rows*cols*sizeof (double);
      b = CreateFullMatrix (rows, cols);
   }   

   memcpy ((void *) &(b -> data[1][1]),(void *) &(a -> data[1][1]), size);

   return b;
}

Matrix MakeFullFromCompact (A)
   Matrix	A;
{
   unsigned 	i,j;
   Matrix	B;

   B = CreateFullMatrix (Mrows(A), Mcols(A));

   for (i = 1 ; i <= Mrows(A) ; i++) 
      for (j = 1 ; j <= Mcols(A) ; j++)
         B -> data [i][j] = mdata(A, i, j);

   return B; 
}

Matrix MakeCompactFromFull (A)
   Matrix	A;
{
   unsigned	*diag;
   Matrix	compA;
   unsigned	i,j,k,
		curr_row;
   unsigned	rows,cols;
   unsigned	size;
   unsigned	height;

   if (IsCompact(A))
      return NullMatrix;

   if (!IsSquare(A))
      return NullMatrix;

   if (!IsSymmetricMatrix(A))
      return NullMatrix;

   rows = Mrows (A);
   cols = Mcols (A);

   size = 0;
   for (i = 1 ; i <= rows ; i++) {
      if (mdata(A,i,i) == 0) {
         size = 1;
         break;
      }
   }

   if (size)
      return NullMatrix;

   diag = (unsigned *) malloc (sizeof(unsigned) * cols);
   if (diag == NULL)
	Fatal ("unable to create compact matrix");

   diag--;
    
	/*
	 * determine the height of the columns and store in diag
	 */

   for (j = 1 ; j <= cols ; j++) {
      diag [j] = 0;
      for (i = 1 ; i <= rows; i++) {

         if (mdata(A,i,j) != 0.0  && i <= j && (j - (i-1)) > diag [j]) 
            diag [j] = j - (i-1);
      }
   }

   size = 0;
   for (i = 1 ; i <= cols ; i++) 
      size += diag [i];
   
   compA = CreateCompactMatrix (rows, cols, size, NULL);

   diag [1] = 1;
   compA -> data [1][1] = mdata (A, 1, 1);
  
   for (i = 2 ; i <= cols ; i++) {
      height = diag [i];
      diag [i] += diag [i-1];
      curr_row = i - height + 1 ;
      for (k = diag [i] - height + 1 ; k <= diag [i] ; k++) 
         compA -> data [k][1] = mdata (A, curr_row++, i);
   }

   compA -> diag = diag;

   return compA;
} 

/**************************************************************************
 * Function:	ConvertRowColumn
 *
 * Parameters:	row		the row of what would be the full matrix
 *		col		the column of what would be the full matrix
 *		a		matrix with array of diagonal address
 *
 * Return value: a valid address if there is one, 0 if there's no
 *		 need to worry about what would be in that spot of
 *		 the full matrix.
 *
 * Description: given a compact column storage scheme described by diag
 *		and convert row and column into an address into
 *		a compact column matrix representation
 *
 **************************************************************************/

int ConvertRowColumn (row, col, a)
    unsigned	row, col;
    Matrix	a;
{
    unsigned	blanks, address;
    unsigned	height;
    unsigned	temp;
    unsigned	*diag;

    diag = a -> diag;

    if (row > col) {
       temp = col;
       col = row;
       row = temp;
    }

    if (col == 1)
       height = 1;
    else
       height = diag [col] - diag [col - 1];

    blanks = col - height;
    if (row > blanks) 
        address = diag [col] + row - col;
    else
        address = 0;

    return address;
}