matrix.c

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%                  M   M   AAA   TTTTT  RRRR   IIIII  X   X                   %
%                  MM MM  A   A    T    R   R    I     X X                    %
%                  M M M  AAAAA    T    RRRR     I      X                     %
%                  M   M  A   A    T    R R      I     X X                    %
%                  M   M  A   A    T    R  R   IIIII  X   X                   %
%                                                                             %
%                                                                             %
%                         MagickCore Matrix Methods                           %
%                                                                             %
%                            Software Design                                  %
%                                 Cristy                                      %
%                              August 2007                                    %
%                                                                             %
%                                                                             %
%  Copyright @ 1999 ImageMagick Studio LLC, a non-profit organization         %
%  dedicated to making software imaging solutions freely available.           %
%                                                                             %
%  You may not use this file except in compliance with the License.  You may  %
%  obtain a copy of the License at                                            %
%                                                                             %
%    https://imagemagick.org/license/                                         %
%                                                                             %
%  Unless required by applicable law or agreed to in writing, software        %
%  distributed under the License is distributed on an "AS IS" BASIS,          %
%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
%  See the License for the specific language governing permissions and        %
%  limitations under the License.                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
*/
␌
/*
  Include declarations.
*/
#include "MagickCore/studio.h"
#include "MagickCore/blob.h"
#include "MagickCore/blob-private.h"
#include "MagickCore/cache.h"
#include "MagickCore/exception.h"
#include "MagickCore/exception-private.h"
#include "MagickCore/image-private.h"
#include "MagickCore/matrix.h"
#include "MagickCore/matrix-private.h"
#include "MagickCore/memory_.h"
#include "MagickCore/memory-private.h"
#include "MagickCore/nt-base-private.h"
#include "MagickCore/pixel-accessor.h"
#include "MagickCore/resource_.h"
#include "MagickCore/semaphore.h"
#include "MagickCore/thread-private.h"
#include "MagickCore/utility.h"
#include "MagickCore/utility-private.h"
␌
/*
  Typedef declaration.
*/
struct _MatrixInfo
{
  CacheType
    type;
 
  size_t
    columns,
    rows,
    stride;
 
  MagickSizeType
    length;
 
  MagickBooleanType
    mapped,
    synchronize;
 
  char
    path[MagickPathExtent];
 
  int
    file;
 
  void
    *elements;
 
  SemaphoreInfo
    *semaphore;
 
  size_t
    signature;
};
␌
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   A c q u i r e M a t r i x I n f o                                         %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AcquireMatrixInfo() allocates the ImageInfo structure.
%
%  The format of the AcquireMatrixInfo method is:
%
%      MatrixInfo *AcquireMatrixInfo(const size_t columns,const size_t rows,
%        const size_t stride,ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o columns: the matrix columns.
%
%    o rows: the matrix rows.
%
%    o stride: the matrix stride.
%
%    o exception: return any errors or warnings in this structure.
%
*/
 
#if defined(SIGBUS)
static void MatrixSignalHandler(int magick_unused(status))
{
  magick_unreferenced(status);
  ThrowFatalException(CacheFatalError,"UnableToExtendMatrixCache");
}
#endif
 
static inline MagickOffsetType WriteMatrixElements(
  const MatrixInfo *magick_restrict matrix_info,const MagickOffsetType offset,
  const MagickSizeType length,const unsigned char *magick_restrict buffer)
{
  MagickOffsetType
    i;
 
  ssize_t
    count;
 
  LockSemaphoreInfo(matrix_info->semaphore);
  if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
    {
      UnlockSemaphoreInfo(matrix_info->semaphore);
      return((MagickOffsetType) -1);
    }
  count=0;
  for (i=0; i < (MagickOffsetType) length; i+=count)
  {
    count=MagickWrite(matrix_info->file,buffer+i,(size_t) MagickMin(length-
      (MagickSizeType) i,(MagickSizeType) MagickMaxBufferExtent));
    if (count <= 0)
      break;
  }
  UnlockSemaphoreInfo(matrix_info->semaphore);
  return(i);
}
 
static MagickBooleanType SetMatrixExtent(
  MatrixInfo *magick_restrict matrix_info,MagickSizeType length)
{
  MagickOffsetType
    count,
    extent,
    offset;
 
  if (length != (MagickSizeType) ((MagickOffsetType) length))
    return(MagickFalse);
  offset=(MagickOffsetType) lseek(matrix_info->file,0,SEEK_END);
  if (offset < 0)
    return(MagickFalse);
  if ((MagickSizeType) offset >= length)
    return(MagickTrue);
  extent=(MagickOffsetType) length-1;
  count=WriteMatrixElements(matrix_info,extent,1,(const unsigned char *) "");
#if defined(MAGICKCORE_HAVE_POSIX_FALLOCATE)
  if (matrix_info->synchronize != MagickFalse)
    (void) posix_fallocate(matrix_info->file,offset+1,extent-offset);
#endif
#if defined(SIGBUS)
  (void) signal(SIGBUS,MatrixSignalHandler);
#endif
  return(count != (MagickOffsetType) 1 ? MagickFalse : MagickTrue);
}
 
MagickExport MatrixInfo *AcquireMatrixInfo(const size_t columns,
  const size_t rows,const size_t stride,ExceptionInfo *exception)
{
  char
    *synchronize;
 
  MagickBooleanType
    status;
 
  MatrixInfo
    *matrix_info;
 
  matrix_info=(MatrixInfo *) AcquireMagickMemory(sizeof(*matrix_info));
  if (matrix_info == (MatrixInfo *) NULL)
    return((MatrixInfo *) NULL);
  (void) memset(matrix_info,0,sizeof(*matrix_info));
  matrix_info->signature=MagickCoreSignature;
  matrix_info->columns=columns;
  matrix_info->rows=rows;
  matrix_info->stride=stride;
  matrix_info->semaphore=AcquireSemaphoreInfo();
  synchronize=GetEnvironmentValue("MAGICK_SYNCHRONIZE");
  if (synchronize != (const char *) NULL)
    {
      matrix_info->synchronize=IsStringTrue(synchronize);
      synchronize=DestroyString(synchronize);
    }
  matrix_info->length=(MagickSizeType) columns*rows*stride;
  if (matrix_info->columns != (size_t) (matrix_info->length/rows/stride))
    {
      (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
        "CacheResourcesExhausted","`%s'","matrix cache");
      return(DestroyMatrixInfo(matrix_info));
    }
  matrix_info->type=MemoryCache;
  status=AcquireMagickResource(AreaResource,matrix_info->length);
  if ((status != MagickFalse) &&
      (matrix_info->length == (MagickSizeType) ((size_t) matrix_info->length)) &&
      ((size_t) matrix_info->length <= GetMaxMemoryRequest()))
    {
      status=AcquireMagickResource(MemoryResource,matrix_info->length);
      if (status != MagickFalse)
        {
          matrix_info->mapped=MagickFalse;
          matrix_info->elements=MagickAssumeAligned(AcquireAlignedMemory(1,
            (size_t) matrix_info->length));
          if (matrix_info->elements == NULL)
            {
              matrix_info->mapped=MagickTrue;
              matrix_info->elements=MapBlob(-1,IOMode,0,(size_t)
                matrix_info->length);
            }
          if (matrix_info->elements == (unsigned short *) NULL)
            RelinquishMagickResource(MemoryResource,matrix_info->length);
        }
    }
  matrix_info->file=(-1);
  if (matrix_info->elements == (unsigned short *) NULL)
    {
      status=AcquireMagickResource(DiskResource,matrix_info->length);
      if (status == MagickFalse)
        {
          (void) ThrowMagickException(exception,GetMagickModule(),CacheError,
            "CacheResourcesExhausted","`%s'","matrix cache");
          return(DestroyMatrixInfo(matrix_info));
        }
      matrix_info->type=DiskCache;
      matrix_info->file=AcquireUniqueFileResource(matrix_info->path);
      if (matrix_info->file == -1)
        return(DestroyMatrixInfo(matrix_info));
      status=AcquireMagickResource(MapResource,matrix_info->length);
      if (status != MagickFalse)
        {
          status=SetMatrixExtent(matrix_info,matrix_info->length);
          if (status != MagickFalse)
            matrix_info->elements=(void *) MapBlob(matrix_info->file,IOMode,0,
              (size_t) matrix_info->length);
          if (matrix_info->elements != NULL)
            matrix_info->type=MapCache;
          else
            RelinquishMagickResource(MapResource,matrix_info->length);
        }
    }
  return(matrix_info);
}
␌
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   A c q u i r e M a g i c k M a t r i x                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AcquireMagickMatrix() allocates and returns a matrix in the form of an
%  array of pointers to an array of doubles, with all values pre-set to zero.
%
%  This used to generate the two dimensional matrix, and vectors required
%  for the GaussJordanElimination() method below, solving some system of
%  simultaneous equations.
%
%  The format of the AcquireMagickMatrix method is:
%
%      double **AcquireMagickMatrix(const size_t number_rows,
%        const size_t size)
%
%  A description of each parameter follows:
%
%    o number_rows: the number pointers for the array of pointers
%      (first dimension).
%
%    o size: the size of the array of doubles each pointer points to
%      (second dimension).
%
*/
MagickExport double **AcquireMagickMatrix(const size_t number_rows,
  const size_t size)
{
  double
    **matrix;
 
  ssize_t
    i,
    j;
 
  matrix=(double **) AcquireQuantumMemory(number_rows,sizeof(*matrix));
  if (matrix == (double **) NULL)
    return((double **) NULL);
  for (i=0; i < (ssize_t) number_rows; i++)
  {
    matrix[i]=(double *) AcquireQuantumMemory(size,sizeof(*matrix[i]));
    if (matrix[i] == (double *) NULL)
      {
        for (j=0; j < i; j++)
          matrix[j]=(double *) RelinquishMagickMemory(matrix[j]);
        matrix=(double **) RelinquishMagickMemory(matrix);
        return((double **) NULL);
      }
    for (j=0; j < (ssize_t) size; j++)
      matrix[i][j]=0.0;
  }
  return(matrix);
}
␌
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   D e s t r o y M a t r i x I n f o                                         %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  DestroyMatrixInfo() dereferences a matrix, deallocating memory associated
%  with the matrix.
%
%  The format of the DestroyImage method is:
%
%      MatrixInfo *DestroyMatrixInfo(MatrixInfo *matrix_info)
%
%  A description of each parameter follows:
%
%    o matrix_info: the matrix.
%
*/
MagickExport MatrixInfo *DestroyMatrixInfo(MatrixInfo *matrix_info)
{
  assert(matrix_info != (MatrixInfo *) NULL);
  assert(matrix_info->signature == MagickCoreSignature);
  LockSemaphoreInfo(matrix_info->semaphore);
  switch (matrix_info->type)
  {
    case MemoryCache:
    {
      if (matrix_info->mapped == MagickFalse)
        matrix_info->elements=RelinquishAlignedMemory(matrix_info->elements);
      else
        {
          (void) UnmapBlob(matrix_info->elements,(size_t) matrix_info->length);
          matrix_info->elements=(unsigned short *) NULL;
        }
      RelinquishMagickResource(MemoryResource,matrix_info->length);
      break;
    }
    case MapCache:
    {
      (void) UnmapBlob(matrix_info->elements,(size_t) matrix_info->length);
      matrix_info->elements=NULL;
      RelinquishMagickResource(MapResource,matrix_info->length);
      magick_fallthrough;
    }
    case DiskCache:
    {
      if (matrix_info->file != -1)
        (void) close_utf8(matrix_info->file);
      (void) RelinquishUniqueFileResource(matrix_info->path);
      RelinquishMagickResource(DiskResource,matrix_info->length);
      break;
    }
    default:
      break;
  }
  UnlockSemaphoreInfo(matrix_info->semaphore);
  RelinquishSemaphoreInfo(&matrix_info->semaphore);
  return((MatrixInfo *) RelinquishMagickMemory(matrix_info));
}
␌
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   G e t M a t r i x C o l u m n s                                           %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  GetMatrixColumns() returns the number of columns in the matrix.
%
%  The format of the GetMatrixColumns method is:
%
%      size_t GetMatrixColumns(const MatrixInfo *matrix_info)
%
%  A description of each parameter follows:
%
%    o matrix_info: the matrix.
%
*/
MagickExport size_t GetMatrixColumns(const MatrixInfo *matrix_info)
{
  assert(matrix_info != (MatrixInfo *) NULL);
  assert(matrix_info->signature == MagickCoreSignature);
  return(matrix_info->columns);
}
␌
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   G e t M a t r i x E l e m e n t                                           %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  GetMatrixElement() returns the specified element in the matrix.
%
%  The format of the GetMatrixElement method is:
%
%      MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info,
%        const ssize_t x,const ssize_t y,void *value)
%
%  A description of each parameter follows:
%
%    o matrix_info: the matrix columns.
%
%    o x: the matrix x-offset.
%
%    o y: the matrix y-offset.
%
%    o value: return the matrix element in this buffer.
%
*/
 
static inline ssize_t EdgeX(const ssize_t x,const size_t columns)
{
  if (x < 0L)
    return(0L);
  if (x >= (ssize_t) columns)
    return((ssize_t) (columns-1));
  return(x);
}
 
static inline ssize_t EdgeY(const ssize_t y,const size_t rows)
{
  if (y < 0L)
    return(0L);
  if (y >= (ssize_t) rows)
    return((ssize_t) (rows-1));
  return(y);
}
 
static inline MagickOffsetType ReadMatrixElements(
  const MatrixInfo *magick_restrict matrix_info,const MagickOffsetType offset,
  const MagickSizeType length,unsigned char *magick_restrict buffer)
{
  MagickOffsetType
    i;
 
  ssize_t
    count;
 
  LockSemaphoreInfo(matrix_info->semaphore);
  if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
    {
      UnlockSemaphoreInfo(matrix_info->semaphore);
      return((MagickOffsetType) -1);
    }
  count=0;
  for (i=0; i < (MagickOffsetType) length; i+=count)
  {
    count=MagickRead(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
      (MagickSizeType) MagickMaxBufferExtent));
    if (count <= 0)
      break;
  }
  UnlockSemaphoreInfo(matrix_info->semaphore);
  return(i);
}
 
MagickExport MagickBooleanType GetMatrixElement(const MatrixInfo *matrix_info,
  const ssize_t x,const ssize_t y,void *value)
{
  MagickOffsetType
    count,
    i;
 
  assert(matrix_info != (const MatrixInfo *) NULL);
  assert(matrix_info->signature == MagickCoreSignature);
  i=EdgeY(y,matrix_info->rows)*(MagickOffsetType) matrix_info->columns+
    EdgeX(x,matrix_info->columns);
  if (matrix_info->type != DiskCache)
    {
      (void) memcpy(value,(unsigned char *) matrix_info->elements+i*
        (MagickOffsetType) matrix_info->stride,matrix_info->stride);
      return(MagickTrue);
    }
  count=ReadMatrixElements(matrix_info,i*(MagickOffsetType) matrix_info->stride,
    matrix_info->stride,(unsigned char *) value);
  if (count != (MagickOffsetType) matrix_info->stride)
    return(MagickFalse);
  return(MagickTrue);
}
␌
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   G e t M a t r i x R o w s                                                 %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  GetMatrixRows() returns the number of rows in the matrix.
%
%  The format of the GetMatrixRows method is:
%
%      size_t GetMatrixRows(const MatrixInfo *matrix_info)
%
%  A description of each parameter follows:
%
%    o matrix_info: the matrix.
%
*/
MagickExport size_t GetMatrixRows(const MatrixInfo *matrix_info)
{
  assert(matrix_info != (const MatrixInfo *) NULL);
  assert(matrix_info->signature == MagickCoreSignature);
  return(matrix_info->rows);
}
␌
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   L e a s t S q u a r e s A d d T e r m s                                   %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  LeastSquaresAddTerms() adds one set of terms and associate results to the
%  given matrix and vectors for solving using least-squares function fitting.
%
%  The format of the AcquireMagickMatrix method is:
%
%      void LeastSquaresAddTerms(double **matrix,double **vectors,
%        const double *terms,const double *results,const size_t rank,
%        const size_t number_vectors);
%
%  A description of each parameter follows:
%
%    o matrix: the square matrix to add given terms/results to.
%
%    o vectors: the result vectors to add terms/results to.
%
%    o terms: the pre-calculated terms (without the unknown coefficient
%             weights) that forms the equation being added.
%
%    o results: the result(s) that should be generated from the given terms
%               weighted by the yet-to-be-solved coefficients.
%
%    o rank: the rank or size of the dimensions of the square matrix.
%            Also the length of vectors, and number of terms being added.
%
%    o number_vectors: Number of result vectors, and number or results being
%      added.  Also represents the number of separable systems of equations
%      that is being solved.
%
%  Example of use...
%
%     2 dimensional Affine Equations (which are separable)
%         c0*x + c2*y + c4*1 => u
%         c1*x + c3*y + c5*1 => v
%
%     double **matrix = AcquireMagickMatrix(3UL,3UL);
%     double **vectors = AcquireMagickMatrix(2UL,3UL);
%     double terms[3], results[2];
%     ...
%     for each given x,y -> u,v
%        terms[0] = x;
%        terms[1] = y;
%        terms[2] = 1;
%        results[0] = u;
%        results[1] = v;
%        LeastSquaresAddTerms(matrix,vectors,terms,results,3UL,2UL);
%     ...
%     if ( GaussJordanElimination(matrix,vectors,3UL,2UL) ) {
%       c0 = vectors[0][0];
%       c2 = vectors[0][1];
%       c4 = vectors[0][2];
%       c1 = vectors[1][0];
%       c3 = vectors[1][1];
%       c5 = vectors[1][2];
%     }
%     else
%       printf("Matrix unsolvable\n");
%     RelinquishMagickMatrix(matrix,3UL);
%     RelinquishMagickMatrix(vectors,2UL);
%
*/
MagickPrivate void LeastSquaresAddTerms(double **matrix,double **vectors,
  const double *terms,const double *results,const size_t rank,
  const size_t number_vectors)
{
  ssize_t
    i,
    j;
 
  for (j=0; j < (ssize_t) rank; j++)
  {
    for (i=0; i < (ssize_t) rank; i++)
      matrix[i][j]+=terms[i]*terms[j];
    for (i=0; i < (ssize_t) number_vectors; i++)
      vectors[i][j]+=results[i]*terms[j];
  }
}
␌
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   G a u s s J o r d a n E l i m n a t i o n                                 %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  GaussJordanElimination() returns a matrix in reduced row echelon form,
%  while simultaneously reducing and thus solving the augmented results
%  matrix.
%
%  The format of the GaussJordanElimination method is:
%
%      MagickBooleanType GaussJordanElimination(double **matrix,
%        double **vectors,const size_t rank,const size_t number_vectors)
%
%  A description of each parameter follows:
%
%    o matrix: the matrix to be reduced, as an 'array of row pointers'.
%
%    o vectors: the additional matrix argumenting the matrix for row reduction.
%             Producing an 'array of column vectors'.
%
%    o rank:  The size of the matrix (both rows and columns).
%             Also represents the number terms that need to be solved.
%
%    o number_vectors: Number of vectors columns, argumenting the above matrix.
%             Usually 1, but can be more for more complex equation solving.
%
%  Note that the 'matrix' is given as a 'array of row pointers' of rank size.
%  That is values can be assigned as   matrix[row][column]   where 'row' is
%  typically the equation, and 'column' is the term of the equation.
%  That is the matrix is in the form of a 'row first array'.
%
%  However 'vectors' is a 'array of column pointers' which can have any number
%  of columns, with each column array the same 'rank' size as 'matrix'.
%
%  This allows for simpler handling of the results, especially is only one
%  column 'vector' is all that is required to produce the desired solution.
%
%  For example, the 'vectors' can consist of a pointer to a simple array of
%  doubles.  when only one set of simultaneous equations is to be solved from
%  the given set of coefficient weighted terms.
%
%     double **matrix = AcquireMagickMatrix(8UL,8UL);
%     double coefficients[8];
%     ...
%     GaussJordanElimination(matrix,&coefficients,8UL,1UL);
%
%  However by specifying more 'columns' (as an 'array of vector columns', you
%  can use this function to solve a set of 'separable' equations.
%
%  For example a distortion function where u = U(x,y) v = V(x,y)
%  And the functions U() and V() have separate coefficients, but are being
%  generated from a common x,y->u,v  data set.
%
%  Another example is generation of a color gradient from a set of colors at
%  specific coordinates, such as a list x,y -> r,g,b,a.
%
%  You can also use the 'vectors' to generate an inverse of the given 'matrix'
%  though as a 'column first array' rather than a 'row first array'. For
%  details see http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
%
*/
MagickPrivate MagickBooleanType GaussJordanElimination(double **matrix,
  double **vectors,const size_t rank,const size_t number_vectors)
{
#define GaussJordanSwap(x,y) \
{ \
  double temp = (x); \
  (x)=(y); \
  (y)=temp; \
}
#define GaussJordanSwapLD(x,y) \
{ \
  long double temp = (x); \
  (x)=(y); \
  (y)=temp; \
}
#define ThrowGaussJordanException() \
{ \
  for (i=0; i < (ssize_t) rank; i++) \
    hp_matrix[i]=(long double *) RelinquishMagickMemory(hp_matrix[i]); \
  hp_matrix=(long double **) RelinquishMagickMemory(hp_matrix); \
  if (pivots != (ssize_t *) NULL) \
    pivots=(ssize_t *) RelinquishMagickMemory(pivots); \
  if (rows != (ssize_t *) NULL) \
    rows=(ssize_t *) RelinquishMagickMemory(rows); \
  if (columns != (ssize_t *) NULL) \
    columns=(ssize_t *) RelinquishMagickMemory(columns); \
  return(MagickFalse); \
}
 
  long double
    **hp_matrix = (long double **) NULL,
    scale;
 
  ssize_t
    column,
    *columns = (ssize_t *) NULL,
    i,
    j,
    *pivots = (ssize_t *) NULL,
    row,
    *rows = (ssize_t *) NULL;
 
  /*
    Allocate high precision matrix.
  */
  hp_matrix=(long double **) AcquireQuantumMemory(rank,sizeof(*hp_matrix));
  if (hp_matrix == (long double **) NULL)
    return(MagickFalse);
  for (i=0; i < (ssize_t) rank; i++)
  {
    hp_matrix[i]=(long double *) AcquireQuantumMemory(rank,
      sizeof(*hp_matrix[i]));
    if (hp_matrix[i] == (long double *) NULL)
      ThrowGaussJordanException();
    for (j=0; j < (ssize_t) rank; j++)
      hp_matrix[i][j]=(long double)matrix[i][j];
  }
  columns=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*columns));
  rows=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*rows));
  pivots=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*pivots));
  if ((columns == (ssize_t *) NULL) || (rows == (ssize_t *) NULL) ||
      (pivots == (ssize_t *) NULL))
    ThrowGaussJordanException();
  (void) memset(columns,0,rank*sizeof(*columns));
  (void) memset(rows,0,rank*sizeof(*rows));
  (void) memset(pivots,0,rank*sizeof(*pivots));
  for (i=0; i < (ssize_t) rank; i++)
  {
    long double
      max = 0.0;
 
    ssize_t
      k;
 
    /*
      Partial pivoting: find the largest absolute value in the unreduced
      submatrix.
    */
    column=(-1);
    row=(-1);
    for (j=0; j < (ssize_t) rank; j++)
      if (pivots[j] != 1)
        for (k=0; k < (ssize_t) rank; k++)
          if ((pivots[k] == 0) && (fabsl(hp_matrix[j][k]) > max))
            {
              max=fabsl(hp_matrix[j][k]);
              row=j;
              column=k;
            }
    if ((column == -1) || (row == -1) || (fabsl(max) < LDBL_MIN))
      ThrowGaussJordanException();  /* Singular or nearly singular matrix */
    pivots[column]++;
    if (row != column)
      {
        for (k=0; k < (ssize_t) rank; k++)
          GaussJordanSwapLD(hp_matrix[row][k],hp_matrix[column][k]);
        for (k=0; k < (ssize_t) number_vectors; k++)
          GaussJordanSwap(vectors[k][row],vectors[k][column]);
      }
    rows[i]=row;
    columns[i]=column;
    if (fabsl(hp_matrix[column][column]) < LDBL_MIN)
      ThrowGaussJordanException();  /* Singular matrix */
    scale=1.0L/hp_matrix[column][column];
    hp_matrix[column][column]=1.0;
    for (j=0; j < (ssize_t) rank; j++)
      hp_matrix[column][j]*=scale;
    for (j=0; j < (ssize_t) number_vectors; j++)
      vectors[j][column]*=(double) scale;
    for (j=0; j < (ssize_t) rank; j++)
      if (j != column)
        {
          scale=hp_matrix[j][column];
          hp_matrix[j][column]=0.0;
          for (k=0; k < (ssize_t) rank; k++)
            hp_matrix[j][k]-=scale*hp_matrix[column][k];
          for (k=0; k < (ssize_t) number_vectors; k++)
            vectors[k][j]-=(double)(scale*(long double) vectors[k][column]);
        }
  }
  for (j=(ssize_t) rank-1; j >= 0; j--)
    if (columns[j] != rows[j])
      for (i=0; i < (ssize_t) rank; i++)
        GaussJordanSwapLD(hp_matrix[i][columns[j]],hp_matrix[i][rows[j]]);
  /*
    Copy back the result to the original matrix.
  */
  for (i=0; i < (ssize_t) rank; i++)
    for (j=0; j < (ssize_t) rank; j++)
      matrix[i][j]=(double)hp_matrix[i][j];
  /*
    Free resources.
  */
  for (i=0; i < (ssize_t) rank; i++)
    hp_matrix[i]=(long double *) RelinquishMagickMemory(hp_matrix[i]);
  hp_matrix=(long double **) RelinquishMagickMemory(hp_matrix);
  pivots=(ssize_t *) RelinquishMagickMemory(pivots);
  rows=(ssize_t *) RelinquishMagickMemory(rows);
  columns=(ssize_t *) RelinquishMagickMemory(columns);
  return(MagickTrue);
}
␌
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   M a t r i x T o I m a g e                                                 %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  MatrixToImage() returns a matrix as an image.  The matrix elements must be
%  of type double otherwise nonsense is returned.
%
%  The format of the MatrixToImage method is:
%
%      Image *MatrixToImage(const MatrixInfo *matrix_info,
%        ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o matrix_info: the matrix.
%
%    o exception: return any errors or warnings in this structure.
%
*/
MagickExport Image *MatrixToImage(const MatrixInfo *matrix_info,
  ExceptionInfo *exception)
{
  CacheView
    *image_view;
 
  double
    max_value,
    min_value,
    scale_factor;
 
  Image
    *image;
 
  MagickBooleanType
    status;
 
  ssize_t
    y;
 
  assert(matrix_info != (const MatrixInfo *) NULL);
  assert(matrix_info->signature == MagickCoreSignature);
  assert(exception != (ExceptionInfo *) NULL);
  assert(exception->signature == MagickCoreSignature);
  if (matrix_info->stride < sizeof(double))
    return((Image *) NULL);
  /*
    Determine range of matrix.
  */
  (void) GetMatrixElement(matrix_info,0,0,&min_value);
  max_value=min_value;
  for (y=0; y < (ssize_t) matrix_info->rows; y++)
  {
    ssize_t
      x;
 
    for (x=0; x < (ssize_t) matrix_info->columns; x++)
    {
      double
        value;
 
      if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
        continue;
      if (value < min_value)
        min_value=value;
      else
        if (value > max_value)
          max_value=value;
    }
  }
  if ((min_value == 0.0) && (max_value == 0.0))
    scale_factor=0;
  else
    if (min_value == max_value)
      {
        scale_factor=(double) QuantumRange/min_value;
        min_value=0;
      }
    else
      scale_factor=(double) QuantumRange/(max_value-min_value);
  /*
    Convert matrix to image.
  */
  image=AcquireImage((ImageInfo *) NULL,exception);
  image->columns=matrix_info->columns;
  image->rows=matrix_info->rows;
  image->colorspace=GRAYColorspace;
  status=MagickTrue;
  image_view=AcquireAuthenticCacheView(image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(image,image,image->rows,2)
#endif
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    double
      value;
 
    Quantum
      *q;
 
    ssize_t
      x;
 
    if (status == MagickFalse)
      continue;
    q=QueueCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
    if (q == (Quantum *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) image->columns; x++)
    {
      if (GetMatrixElement(matrix_info,x,y,&value) == MagickFalse)
        continue;
      value=scale_factor*(value-min_value);
      *q=ClampToQuantum(value);
      q+=(ptrdiff_t) GetPixelChannels(image);
    }
    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
      status=MagickFalse;
  }
  image_view=DestroyCacheView(image_view);
  if (status == MagickFalse)
    image=DestroyImage(image);
  return(image);
}
␌
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   N u l l M a t r i x                                                       %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  NullMatrix() sets all elements of the matrix to zero.
%
%  The format of the memset method is:
%
%      MagickBooleanType *NullMatrix(MatrixInfo *matrix_info)
%
%  A description of each parameter follows:
%
%    o matrix_info: the matrix.
%
*/
MagickExport MagickBooleanType NullMatrix(MatrixInfo *matrix_info)
{
  ssize_t
    x;
 
  ssize_t
    count,
    y;
 
  unsigned char
    value;
 
  assert(matrix_info != (const MatrixInfo *) NULL);
  assert(matrix_info->signature == MagickCoreSignature);
  if (matrix_info->type != DiskCache)
    {
      (void) memset(matrix_info->elements,0,(size_t)
        matrix_info->length);
      return(MagickTrue);
    }
  value=0;
  (void) lseek(matrix_info->file,0,SEEK_SET);
  for (y=0; y < (ssize_t) matrix_info->rows; y++)
  {
    for (x=0; x < (ssize_t) matrix_info->length; x++)
    {
      count=MagickWrite(matrix_info->file,&value,sizeof(value));
      if (count != (ssize_t) sizeof(value))
        break;
    }
    if (x < (ssize_t) matrix_info->length)
      break;
  }
  return(y < (ssize_t) matrix_info->rows ? MagickFalse : MagickTrue);
}
␌
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   R e l i n q u i s h M a g i c k M a t r i x                               %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  RelinquishMagickMatrix() frees the previously acquired matrix (array of
%  pointers to arrays of doubles).
%
%  The format of the RelinquishMagickMatrix method is:
%
%      double **RelinquishMagickMatrix(double **matrix,
%        const size_t number_rows)
%
%  A description of each parameter follows:
%
%    o matrix: the matrix to relinquish
%
%    o number_rows: the first dimension of the acquired matrix (number of
%      pointers)
%
*/
MagickExport double **RelinquishMagickMatrix(double **matrix,
  const size_t number_rows)
{
  ssize_t
    i;
 
  if (matrix == (double **) NULL )
    return(matrix);
  for (i=0; i < (ssize_t) number_rows; i++)
    matrix[i]=(double *) RelinquishMagickMemory(matrix[i]);
  matrix=(double **) RelinquishMagickMemory(matrix);
  return(matrix);
}
␌
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   S e t M a t r i x E l e m e n t                                           %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  SetMatrixElement() sets the specified element in the matrix.
%
%  The format of the SetMatrixElement method is:
%
%      MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info,
%        const ssize_t x,const ssize_t y,void *value)
%
%  A description of each parameter follows:
%
%    o matrix_info: the matrix columns.
%
%    o x: the matrix x-offset.
%
%    o y: the matrix y-offset.
%
%    o value: set the matrix element to this value.
%
*/
 
MagickExport MagickBooleanType SetMatrixElement(const MatrixInfo *matrix_info,
  const ssize_t x,const ssize_t y,const void *value)
{
  MagickOffsetType
    count,
    i;
 
  assert(matrix_info != (const MatrixInfo *) NULL);
  assert(matrix_info->signature == MagickCoreSignature);
  i=y*(MagickOffsetType) matrix_info->columns+x;
  if ((i < 0) ||
      (((MagickSizeType) i*matrix_info->stride) >= matrix_info->length))
    return(MagickFalse);
  if (matrix_info->type != DiskCache)
    {
      (void) memcpy((unsigned char *) matrix_info->elements+i*
        (MagickOffsetType) matrix_info->stride,value,matrix_info->stride);
      return(MagickTrue);
    }
  count=WriteMatrixElements(matrix_info,i*(MagickOffsetType)
    matrix_info->stride,matrix_info->stride,(unsigned char *) value);
  if (count != (MagickOffsetType) matrix_info->stride)
    return(MagickFalse);
  return(MagickTrue);
}