 /*
  * Khoros: $Id: lmmult.c,v 1.3 1992/03/20 23:47:08 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lmmult.c,v 1.3 1992/03/20 23:47:08 dkhoros Exp $";
#endif

 /*
  * $Log: lmmult.c,v $
 * Revision 1.3  1992/03/20  23:47:08  dkhoros
 * VirtualPatch5
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1992, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as to the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including, for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1992 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: lmmult.c
 >>>>
 >>>>      Program Name: mmult
 >>>>
 >>>> Date Last Updated: Thu Mar  5 08:53:10 1992 
 >>>>
 >>>>          Routines: lmmult - the library call for mmult
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lmmult - library call for mmult
*
* Purpose:
*    
*    mmult - Multiplies two matricies
*    
*    

* Input:
*    
*    image1         the first matrix in multiplication sequence
*    
*    image2         the second matrix in the multiplication sequence
*    
*    

* Output:
*    
*    image1         holds the result of image1*image2, the  dimensions
*                   should be rows1*cols2
*    
*    Return Value:  1 on success, 0 on failure.
*    
*    

*
* Written By: Matthew Lawrence, Jeremy Worley
*    
*    Jeremy Worley 05 Mar 1992 08:52 MST
*              Made certain implicitly declared  functions  explicitly
*              declared.
*    
*    

****************************************************************/


/* -library_def */
int lmmult(image1,image2)
   struct xvimage *image1,*image2;
/* -library_def_end */

/* -library_code */
{
  int i,rows1,cols1,rows2,cols2,bands,type;
  double *mat,*m1,*m2;
  char *program = "lmmult";
  struct xvimage *images[2];
  int proper_num_images(), match_num_bands(), lvcast(), lmult();

  images[0] = image1;
  images[1] = image2;
 
  rows1 = (int)image1->col_size;
  cols1 = (int)image1->row_size;
  rows2 = (int)image2->col_size;
  cols2 = (int)image2->row_size;

/*
** check some stuff out
*/

  if(!proper_num_images(program,image1,1,FALSE)){
     fprintf(stderr,"%s:  1 image structures should have only one image.\n",
         program);
     return(0);
  }

  if(!proper_num_images(program,image2,1,FALSE)){
     fprintf(stderr,"%s:  2 image structures should have only one image.\n",
         program);
     return(0);
  }

  if(!match_num_bands(program,image1,image2,FALSE)){
     fprintf(stderr,"%s:  images structures must have same number of bands.\n",
         program);
     return(0);
  }

  bands = image1->num_data_bands;

/*
** cast all of the data up to double or double complex
*/

   if(image1->data_storage_type==VFF_TYP_COMPLEX || image1->data_storage_type==
        VFF_TYP_DCOMPLEX || image2->data_storage_type==VFF_TYP_COMPLEX ||
        image2->data_storage_type==VFF_TYP_DCOMPLEX){
      if(!lvcast(images,(int)2,(int)VFF_TYP_DCOMPLEX)){
         fprintf(stderr,"%s:  Failure to cast up to double complex.\n",program);
         return(0);
      }
      type = 1;
   }else{
      if(!lvcast(images,(int)2,(int)VFF_TYP_DOUBLE)){
         fprintf(stderr,"%s:  Failure to cast up to double.\n",program);
         return(0);
      }
      type = 0;
   }

   if((mat=(double *)malloc((type+1)*bands*cols2*rows1*sizeof(double)))==NULL){
      fprintf(stderr,"%s:  failure to malloc.\n",program);
      return(0);
   }

/*
** process all of the data
*/

   m1 = (double *)(image1->imagedata);
   m2 = (double *)(image2->imagedata);

   for(i=0;i<bands;i++){
      if(!lmult(&mat[i*rows1*cols2],m1,rows1,cols1,type,m2,rows2,cols2,type)){
         fprintf(stderr,"%s:  Failed call to static routine.\n",program);
         free((char *)mat);
         return(0);
      }

      m1 += (type+1)*rows1*cols1;
      m2 += (type+1)*rows2*cols2;
   }

   image1->imagedata = (char *)mat;
   image1->row_size = cols2;
   image1->col_size = rows1;

   return(1);
}
/***********************************************************************
*
*  Routine Name: 
*
*          Date:
*        
*       Purpose:  
*
*         Input: 
*
*        Output: 
*
*    Written By:  
*
* Modifications:
*
***********************************************************************/

int lmult(m3,m1,rows1,cols1,type1,m2,rows2,cols2,type2)
   double *m3,*m1,*m2;
   int rows1,cols1,type1,rows2,cols2,type2;
   {
   char *program = "lmmult";
   double *temp,*store;
   int tempcount,i,temprc;
   int m1pos,m2pos,m3pos,m3postemp,r1,r2,c2;
   if(cols1 != rows2)
      {
      (void)fprintf(stderr,"%s: number of columns of first matrix",program);
      (void)fprintf(stderr," must equal number of rows of second matrix to multiply\n");
      return(0);
      }

   if(rows1==0 || cols1==0 || rows2==0 || cols2==0)
      {
      (void)fprintf(stderr,"%s: zero dimension matrices not allowed\n",
         program);
      return(0);
      }

   if(m3==NULL || m2==NULL || m1==NULL)
      {
      (void)fprintf(stderr,"%s: Null input/output parameter not allowed.\n",program);
      return(0);
      }

   if(type1 < 0 || type1 > 1 || type2 < 0 || type2 > 1)
      {
      (void)fprintf(stderr,"%s: illegal type paramter specified \n",program);
      return(0);
      }

   /* 
    *  Make matricies the same type.
    */
   if(type1 != type2)
      {   /* One matix is complex, convert the other matrix. */
      temp = (double *)malloc(sizeof(double)*(rows2*cols2*2));
      if(temp == NULL)
         {
         (void)fprintf(stderr,"%s: unable to allocate memory.\n",program);
         return(0);
         }
      if(type1 == 1)
         temprc = rows2*cols2;    /* Matrix 2 needs to be converted. */
      else
         temprc = rows1*cols1;    /* Matrix 1 needs to be converted. */
      tempcount = 0;
      for(i = 0; i < temprc; i++)
         {
         if(type1 == 1)
            *(temp+tempcount) = *(m2 + i);
         else
            *(temp+tempcount) = *(m1 + i);
         *(temp+tempcount+1) = 0;
         tempcount = tempcount + 2;
         }
      if(type1 == 1)
         {
         store = m2;
         m2 = temp;
         temp = store;
         free(temp);
         }
      else
         {
         store = m1;
         m1 = temp;
         temp = store;
         free(temp);
         }
      type1 = 1;
      /* type2 = 1; */
      }


   for(r1 = 0; r1 < (type1+1)*rows1*cols2; r1++)  
      *(m3+r1) = 0; /* Initialize the result matrix. */
   m1pos = 0;
   m2pos = 0;
   m3pos = 0;
   /* m3postemp = 0; */

   /*
    * Start of multiply
    */

   for(r1 = 0; r1 < rows1; r1++)
      {
      for(r2 = 0; r2 < rows2; r2++)
         {
         m3postemp = m3pos;
         for(c2 = 0; c2 < cols2; c2++)
            {
            if(type1 == 1)
               {
               *(m3+m3postemp) = (*(m1+m1pos)) * (*(m2+m2pos)) - (*(m1+m1pos+1)) *
                  (*(m2+m2pos+1)) + *(m3+m3postemp); 
               *(m3+m3postemp+1) = (*(m1+m1pos)) * (*(m2+m2pos+1)) + (*(m1+m1pos+1)) *
                  (*(m2+m2pos)) + *(m3+m3postemp+1);
               m2pos = m2pos + 2;
               m3postemp = m3postemp + 2;
               }
            else
               {
               *(m3+m3postemp) = (*(m1+m1pos)) * (*(m2+m2pos)) + *(m3+m3postemp);
               m2pos++;
               m3postemp++;
               }
            } /* End of c2 */
         if(type1 == 1)
            m1pos = m1pos + 2;
         else
            m1pos++;
         } /* End of r2 */
      m2pos = 0;
      if(type1 == 1)
         m3pos = m3pos + cols2 * 2;
      else
         m3pos = m3pos + cols2;
      } /* End of r1 */

      return(1);  /* Return a valid exit value. */
   } /* End of library routine */

/* -library_code_end */
