text: /*
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                                                                             %
  %                                                                             %
  %                                                                             %
  %                  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-2020 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/script/license.php                               %
  %                                                                             %
  %  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/pixel-accessor.h"
  #include "MagickCore/pixel-private.h"
  #include "MagickCore/resource_.h"
  #include "MagickCore/semaphore.h"
  #include "MagickCore/thread-private.h"
  #include "MagickCore/utility.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 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;
  
  #if !defined(MAGICKCORE_HAVE_PWRITE)
    LockSemaphoreInfo(matrix_info->semaphore);
    if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
      {
        UnlockSemaphoreInfo(matrix_info->semaphore);
        return((MagickOffsetType) -1);
      }
  #endif
    count=0;
    for (i=0; i < (MagickOffsetType) length; i+=count)
    {
  #if !defined(MAGICKCORE_HAVE_PWRITE)
      count=write(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
        (MagickSizeType) SSIZE_MAX));
  #else
      count=pwrite(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
        (MagickSizeType) SSIZE_MAX),(off_t) (offset+i));
  #endif
      if (count <= 0)
        {
          count=0;
          if (errno != EINTR)
            break;
        }
    }
  #if !defined(MAGICKCORE_HAVE_PWRITE)
    UnlockSemaphoreInfo(matrix_info->semaphore);
  #endif
    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)))
      {
        status=AcquireMagickResource(MemoryResource,matrix_info->length);
        if (status != MagickFalse)
          {
            matrix_info->mapped=MagickFalse;
            matrix_info->elements=AcquireMagickMemory((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
  %  simultanious 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=RelinquishMagickMemory(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);
      }
      case DiskCache:
      {
        if (matrix_info->file != -1)
          (void) close(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 a u s s J o r d a n E l i m i n a t i o n                               %
  %                                                                             %
  %                                                                             %
  %                                                                             %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  %  GaussJordanElimination() returns a matrix in reduced row echelon form,
  %  while simultaneously reducing and thus solving the augumented results
  %  matrix.
  %
  %  See also  http://en.wikipedia.org/wiki/Gauss-Jordan_elimination
  %
  %  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.
  %             Usally 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 simultanious equations is to be solved from
  %  the given set of coefficient weighted terms.
  %
  %     double **matrix = AcquireMagickMatrix(8UL,8UL);
  %     double coefficents[8];
  %     ...
  %     GaussJordanElimination(matrix, &coefficents, 8UL, 1UL);
  %
  %  However by specifing 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 coefficents, 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 coordients, 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) \
  { \
    if ((x) != (y)) \
      { \
        (x)+=(y); \
        (y)=(x)-(y); \
        (x)=(x)-(y); \
      } \
  }
  
    double
      max,
      scale;
  
    ssize_t
      i,
      j,
      k;
  
    ssize_t
      column,
      *columns,
      *pivots,
      row,
      *rows;
  
    columns=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*columns));
    rows=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*rows));
    pivots=(ssize_t *) AcquireQuantumMemory(rank,sizeof(*pivots));
    if ((rows == (ssize_t *) NULL) || (columns == (ssize_t *) NULL) ||
        (pivots == (ssize_t *) NULL))
      {
        if (pivots != (ssize_t *) NULL)
          pivots=(ssize_t *) RelinquishMagickMemory(pivots);
        if (columns != (ssize_t *) NULL)
          columns=(ssize_t *) RelinquishMagickMemory(columns);
        if (rows != (ssize_t *) NULL)
          rows=(ssize_t *) RelinquishMagickMemory(rows);
        return(MagickFalse);
      }
    (void) memset(columns,0,rank*sizeof(*columns));
    (void) memset(rows,0,rank*sizeof(*rows));
    (void) memset(pivots,0,rank*sizeof(*pivots));
    column=0;
    row=0;
    for (i=0; i < (ssize_t) rank; i++)
    {
      max=0.0;
      for (j=0; j < (ssize_t) rank; j++)
        if (pivots[j] != 1)
          {
            for (k=0; k < (ssize_t) rank; k++)
              if (pivots[k] != 0)
                {
                  if (pivots[k] > 1)
                    return(MagickFalse);
                }
              else
                if (fabs(matrix[j][k]) >= max)
                  {
                    max=fabs(matrix[j][k]);
                    row=j;
                    column=k;
                  }
          }
      pivots[column]++;
      if (row != column)
        {
          for (k=0; k < (ssize_t) rank; k++)
            GaussJordanSwap(matrix[row][k],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 (matrix[column][column] == 0.0)
        return(MagickFalse);  /* sigularity */
      scale=PerceptibleReciprocal(matrix[column][column]);
      matrix[column][column]=1.0;
      for (j=0; j < (ssize_t) rank; j++)
        matrix[column][j]*=scale;
      for (j=0; j < (ssize_t) number_vectors; j++)
        vectors[j][column]*=scale;
      for (j=0; j < (ssize_t) rank; j++)
        if (j != column)
          {
            scale=matrix[j][column];
            matrix[j][column]=0.0;
            for (k=0; k < (ssize_t) rank; k++)
              matrix[j][k]-=scale*matrix[column][k];
            for (k=0; k < (ssize_t) number_vectors; k++)
              vectors[k][j]-=scale*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++)
          GaussJordanSwap(matrix[i][rows[j]],matrix[i][columns[j]]);
    pivots=(ssize_t *) RelinquishMagickMemory(pivots);
    rows=(ssize_t *) RelinquishMagickMemory(rows);
    columns=(ssize_t *) RelinquishMagickMemory(columns);
    return(MagickTrue);
  }
  
  /*
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                                                                             %
  %                                                                             %
  %                                                                             %
  %   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 specifed 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;
  
  #if !defined(MAGICKCORE_HAVE_PREAD)
    LockSemaphoreInfo(matrix_info->semaphore);
    if (lseek(matrix_info->file,offset,SEEK_SET) < 0)
      {
        UnlockSemaphoreInfo(matrix_info->semaphore);
        return((MagickOffsetType) -1);
      }
  #endif
    count=0;
    for (i=0; i < (MagickOffsetType) length; i+=count)
    {
  #if !defined(MAGICKCORE_HAVE_PREAD)
      count=read(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
        (MagickSizeType) SSIZE_MAX));
  #else
      count=pread(matrix_info->file,buffer+i,(size_t) MagickMin(length-i,
        (MagickSizeType) SSIZE_MAX),(off_t) (offset+i));
  #endif
      if (count <= 0)
        {
          count=0;
          if (errno != EINTR)
            break;
        }
    }
  #if !defined(MAGICKCORE_HAVE_PREAD)
    UnlockSemaphoreInfo(matrix_info->semaphore);
  #endif
    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=(MagickOffsetType) EdgeY(y,matrix_info->rows)*matrix_info->columns+
      EdgeX(x,matrix_info->columns);
    if (matrix_info->type != DiskCache)
      {
        (void) memcpy(value,(unsigned char *) matrix_info->elements+i*
          matrix_info->stride,matrix_info->stride);
        return(MagickTrue);
      }
    count=ReadMatrixElements(matrix_info,i*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 coefficent
  %             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 coefficents.
  %
  %    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];
    }
  }
  
  /*
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                                                                             %
  %                                                                             %
  %                                                                             %
  %   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,1)
  #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+=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=write(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 specifed 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=(MagickOffsetType) y*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*
          matrix_info->stride,value,matrix_info->stride);
        return(MagickTrue);
      }
    count=WriteMatrixElements(matrix_info,i*matrix_info->stride,
      matrix_info->stride,(unsigned char *) value);
    if (count != (MagickOffsetType) matrix_info->stride)
      return(MagickFalse);
    return(MagickTrue);
  }
  