lapack_interface.c

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/* This file is part of SCIPSDP - a solving framework for mixed-integer      */
/* semidefinite programs based on SCIP.                                      */
/*                                                                           */
/* Copyright (C) 2011-2013 Discrete Optimization, TU Darmstadt               */
/*                         EDOM, FAU Erlangen-Nürnberg                       */
/*               2014-2021 Discrete Optimization, TU Darmstadt               */
/*                                                                           */
/*                                                                           */
/* This program is free software; you can redistribute it and/or             */
/* modify it under the terms of the GNU Lesser General Public License        */
/* as published by the Free Software Foundation; either version 3            */
/* 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 Lesser General Public License for more details.                       */
/*                                                                           */
/* You should have received a copy of the GNU Lesser General Public License  */
/* along with this program; if not, write to the Free Software               */
/* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.*/
/*                                                                           */
/*                                                                           */
/* Based on SCIP - Solving Constraint Integer Programs                       */
/* Copyright (C) 2002-2021 Zuse Institute Berlin                             */
/* SCIP is distributed under the terms of the SCIP Academic Licence,         */
/* see file COPYING in the SCIP distribution.                                */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/* #define PRINTMATRICES     /\* Should all matrices appearing in best rank-1 approximation heuristic be printed? *\/ */

#include <assert.h>

#include "lapack_interface.h"
#include "config.h"                          /* for F77_FUNC */

#include "scip/def.h"
#include "scip/pub_message.h"                /* for debug and error message */
#include "blockmemshell/memory.h"
#include "scip/type_retcode.h"

/* turn off lint warnings for whole file: */
/*lint --e{788,818}*/


/* we use 64 bit integers as the base type */
typedef long long int LAPACKINTTYPE;


#define BMS_CALL(x)   do                                                                                      \
                      {                                                                                       \
                          if( NULL == (x) )                                                                   \
                          {                                                                                   \
                             SCIPerrorMessage("No memory in function call\n");                                \
                             return SCIP_NOMEMORY;                                                            \
                          }                                                                                   \
                      }                                                                                       \
                      while( FALSE )

#define SCIP_RealTOINT(x) ((LAPACKINTTYPE) (x + 0.5))

/*
 * BLAS/LAPACK Calls
 */

void F77_FUNC(dsyevr, DSYEVR)(char* JOBZ, char* RANGE, char* UPLO,
   LAPACKINTTYPE* N, SCIP_Real* A, LAPACKINTTYPE* LDA,
   SCIP_Real* VL, SCIP_Real* VU,
   LAPACKINTTYPE* IL, LAPACKINTTYPE* IU,
   SCIP_Real* ABSTOL, LAPACKINTTYPE* M, SCIP_Real* W, SCIP_Real* Z,
   LAPACKINTTYPE* LDZ, LAPACKINTTYPE* ISUPPZ, SCIP_Real* WORK,
   LAPACKINTTYPE* LWORK, LAPACKINTTYPE* IWORK, LAPACKINTTYPE* LIWORK,
   LAPACKINTTYPE* INFO);

void F77_FUNC(dsyevx, DSYEVX)(char* JOBZ, char* RANGE, char* UPLO,
   LAPACKINTTYPE* N, SCIP_Real* A, LAPACKINTTYPE* LDA,
   SCIP_Real* VL, SCIP_Real* VU,
   LAPACKINTTYPE* IL, LAPACKINTTYPE* IU,
   SCIP_Real* ABSTOL, LAPACKINTTYPE* M, SCIP_Real* W, SCIP_Real* Z,
   LAPACKINTTYPE* LDZ, SCIP_Real* WORK, LAPACKINTTYPE* LWORK, LAPACKINTTYPE* IWORK,
   LAPACKINTTYPE* IFAIL, LAPACKINTTYPE* INFO);

void F77_FUNC(dgemv, DGEMV)(char* TRANS, LAPACKINTTYPE* M,
   LAPACKINTTYPE* N, SCIP_Real* ALPHA, SCIP_Real* A, LAPACKINTTYPE* LDA,
   SCIP_Real* X, LAPACKINTTYPE* INCX, SCIP_Real* BETA, SCIP_Real* Y, LAPACKINTTYPE* INCY);

void F77_FUNC(dgemm, DGEMM)(char* TRANSA, char* TRANSB, LAPACKINTTYPE* M, LAPACKINTTYPE* N, LAPACKINTTYPE* K, SCIP_Real* ALPHA,
      SCIP_Real* A, LAPACKINTTYPE* LDA, SCIP_Real* B, LAPACKINTTYPE* LDB, SCIP_Real* BETA, SCIP_Real* C, LAPACKINTTYPE* LDC);

/* LAPACK Fortran subroutine DGELSD */
void F77_FUNC(dgelsd, DGELSD)(LAPACKINTTYPE* M, LAPACKINTTYPE* N, LAPACKINTTYPE* NRHS,
      SCIP_Real* A, LAPACKINTTYPE* LDA, SCIP_Real* b, LAPACKINTTYPE* LDB, SCIP_Real* S, SCIP_Real* RCOND, LAPACKINTTYPE* RANK,
      SCIP_Real* WORK, LAPACKINTTYPE* LWORK, LAPACKINTTYPE* IWORK, LAPACKINTTYPE* INFO);

static
int convertToInt(
   long long int         num                 
   )
{
   long long int work;
   int checkval = 1;

   assert(sizeof(work) > sizeof(checkval)); /*lint !e506*/

   /* if we have a little-endian machine (e.g, x86), the sought value is in the bottom part */
   if ( *(int8_t*)&checkval != 0 ) /*lint !e774*/
   {
      /* if the top part is nonzero, we assume that the number is negative */
      if ( *((int8_t*)&num + 4) != 0 ) /*lint !e2662*/
      {
         work = -num;
         return -(*((int*)&work));
      }
      return *((int*)&num);
   }

   /* otherwise we have a big-endian machine (e.g., PowerPC); the sought value is in the top part */
   assert( *(int8_t*)&checkval == 0 );

   /* if the bottom part is nonzero, we assume that the number is negative */
   if ( *(int8_t*)&num != 0 ) /*lint !e774*/
   {
      work = -num;
      return -(*((int*)&work + 4)); /*lint !e2662*/
   }
   return *((int*)&num + 4);
}


/*
 * Functions
 */

SCIP_RETCODE SCIPlapackComputeIthEigenvalue(
   BMS_BUFMEM*           bufmem,             
   SCIP_Bool             geteigenvectors,    
   int                   n,                  
   SCIP_Real*            A,                  
   int                   i,                  
   SCIP_Real*            eigenvalue,         
   SCIP_Real*            eigenvector         
   )
{
   LAPACKINTTYPE* IWORK;
   LAPACKINTTYPE* ISUPPZ;
   LAPACKINTTYPE N;
   LAPACKINTTYPE INFO;
   LAPACKINTTYPE LDA;
   LAPACKINTTYPE WISIZE;
   LAPACKINTTYPE IL;
   LAPACKINTTYPE IU;
   LAPACKINTTYPE M;
   LAPACKINTTYPE LDZ;
   LAPACKINTTYPE LWORK;
   LAPACKINTTYPE LIWORK;
   SCIP_Real* WORK;
   SCIP_Real* WTMP;
   SCIP_Real ABSTOL;
   SCIP_Real WSIZE;
   SCIP_Real VL;
   SCIP_Real VU;
   char JOBZ;
   char RANGE;
   char UPLO;

   assert( bufmem != NULL );
   assert( n > 0 );
   assert( n < INT_MAX );
   assert( A != NULL );
   assert( 0 < i && i <= n );
   assert( eigenvalue != NULL );
   assert( ! geteigenvectors || eigenvector != NULL );

   N = n;
   JOBZ = geteigenvectors ? 'V' : 'N';
   RANGE = 'I';
   UPLO = 'L';
   LDA  = n;
   ABSTOL = 0.0; /* we use abstol = 0, since some lapack return an error otherwise */
   VL = -1e20;
   VU = 1e20;
   IL = i;
   IU = i;
   M = 1;
   LDZ = n;
   INFO = 0LL;

   /* standard LAPACK workspace query, to get the amount of needed memory */
   LWORK = -1LL;
   LIWORK = -1LL;

   /* this computes the internally needed memory and returns this as (the first entries of [the 1x1 arrays]) WSIZE and WISIZE */
   F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
      &N, NULL, &LDA,
      NULL, NULL,
      &IL, &IU,
      &ABSTOL, &M, NULL, NULL,
      &LDZ, NULL, &WSIZE,
      &LWORK, &WISIZE, &LIWORK,
      &INFO);

   if ( convertToInt(INFO) != 0 )
   {
      SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d.\n", convertToInt(INFO));
      return SCIP_ERROR;
   }

   /* allocate workspace */
   LWORK = SCIP_RealTOINT(WSIZE);
   LIWORK = WISIZE;

   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WORK, (int) LWORK) );
   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &IWORK, (int) LIWORK) );
   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WTMP, (int) N) );
   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &ISUPPZ, 2) ); /*lint !e506*/

   /* call the function */
   F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
      &N, A, &LDA,
      &VL, &VU,
      &IL, &IU,
      &ABSTOL, &M, WTMP, eigenvector,
      &LDZ, ISUPPZ, WORK,
      &LWORK, IWORK, &LIWORK,
      &INFO);

   /* handle output */
   if ( convertToInt(INFO) == 0 )
     *eigenvalue = WTMP[0];

   /* free memory */
   BMSfreeBufferMemoryArray(bufmem, &ISUPPZ);
   BMSfreeBufferMemoryArray(bufmem, &WTMP);
   BMSfreeBufferMemoryArray(bufmem, &IWORK);
   BMSfreeBufferMemoryArray(bufmem, &WORK);

   if ( convertToInt(INFO) != 0 )
   {
      SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d.\n", convertToInt(INFO));
      return SCIP_ERROR;
   }

   return SCIP_OKAY;
}

SCIP_RETCODE SCIPlapackComputeIthEigenvalueAlternative(
   BMS_BUFMEM*           bufmem,             
   SCIP_Bool             geteigenvectors,    
   int                   n,                  
   SCIP_Real*            A,                  
   int                   i,                  
   SCIP_Real*            eigenvalue,         
   SCIP_Real*            eigenvector         
   )
{
   LAPACKINTTYPE* IWORK;
   LAPACKINTTYPE* IFAIL;
   LAPACKINTTYPE N;
   LAPACKINTTYPE INFO;
   LAPACKINTTYPE LDA;
   LAPACKINTTYPE IL;
   LAPACKINTTYPE IU;
   LAPACKINTTYPE M;
   LAPACKINTTYPE LDZ;
   LAPACKINTTYPE LWORK;
   SCIP_Real* WORK;
   SCIP_Real* WTMP;
   SCIP_Real WSIZE;
   SCIP_Real ABSTOL;
   SCIP_Real VL;
   SCIP_Real VU;
   char JOBZ;
   char RANGE;
   char UPLO;

   assert( bufmem != NULL );
   assert( n > 0 );
   assert( n < INT_MAX );
   assert( A != NULL );
   assert( 0 < i && i <= n );
   assert( eigenvalue != NULL );
   assert( ! geteigenvectors || eigenvector != NULL );

   N = n;
   JOBZ = geteigenvectors ? 'V' : 'N';
   RANGE = 'I';
   UPLO = 'L';
   LDA  = n;
   ABSTOL = 0.0;  /* we use abstol = 0, since some lapack return an error otherwise */
   VL = -1e20;
   VU = 1e20;
   IL = i;
   IU = i;
   M = 1;
   LDZ = n;
   INFO = 0LL;

   /* standard LAPACK workspace query, to get the amount of needed memory */
   LWORK = -1LL;

   /* this computes the internally needed memory and returns this as (the first entries of [the 1x1 arrays]) WSIZE */
   F77_FUNC(dsyevx, DSYEVX)( &JOBZ, &RANGE, &UPLO,
      &N, NULL, &LDA,
      NULL, NULL,
      &IL, &IU,
      &ABSTOL, &M, NULL, NULL,
      &LDZ, &WSIZE, &LWORK, NULL, NULL,
      &INFO);

   if ( convertToInt(INFO) != 0 )
   {
      SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d.\n", convertToInt(INFO));
      return SCIP_ERROR;
   }

   /* allocate workspace */
   LWORK = SCIP_RealTOINT(WSIZE);

   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WORK, (int) LWORK) );
   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &IWORK, (int) 5 * N) );
   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WTMP, (int) N) );
   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &IFAIL, (int) N) ); /*lint !e506*/

   /* call the function */
   F77_FUNC(dsyevx, DSYEVX)( &JOBZ, &RANGE, &UPLO,
      &N, A, &LDA,
      &VL, &VU,
      &IL, &IU,
      &ABSTOL, &M, WTMP, eigenvector,
      &LDZ, WORK, &LWORK, IWORK, IFAIL,
      &INFO);

   /* handle output */
   if ( convertToInt(INFO) == 0 )
     *eigenvalue = WTMP[0];

   /* free memory */
   BMSfreeBufferMemoryArray(bufmem, &IFAIL);
   BMSfreeBufferMemoryArray(bufmem, &WTMP);
   BMSfreeBufferMemoryArray(bufmem, &IWORK);
   BMSfreeBufferMemoryArray(bufmem, &WORK);

   if ( convertToInt(INFO) != 0 )
   {
      SCIPerrorMessage("There was an error when calling DSYEVX. INFO = %d.\n", convertToInt(INFO));
      return SCIP_ERROR;
   }

   return SCIP_OKAY;
}

SCIP_RETCODE SCIPlapackComputeEigenvectorsNegative(
   BMS_BUFMEM*           bufmem,             
   int                   n,                  
   SCIP_Real*            A,                  
   SCIP_Real             tol,                
   int*                  neigenvalues,       
   SCIP_Real*            eigenvalues,        
   SCIP_Real*            eigenvectors        
   )
{
   LAPACKINTTYPE* ISUPPZ;
   LAPACKINTTYPE* IWORK;
   LAPACKINTTYPE N;
   LAPACKINTTYPE INFO;
   LAPACKINTTYPE LDA;
   LAPACKINTTYPE LWORK;
   LAPACKINTTYPE LIWORK;
   LAPACKINTTYPE M;
   LAPACKINTTYPE LDZ;
   LAPACKINTTYPE WISIZE;
   SCIP_Real* WORK;
   SCIP_Real ABSTOL;
   SCIP_Real WSIZE;
   SCIP_Real VL;
   SCIP_Real VU;
   char JOBZ;
   char RANGE;
   char UPLO;

   assert( bufmem != NULL );
   assert( n > 0 );
   assert( n < INT_MAX );
   assert( A != NULL );
   assert( tol >= 0 );
   assert( neigenvalues != NULL );
   assert( eigenvalues != NULL );
   assert( eigenvectors != NULL );

   N = n;
   JOBZ = 'V';
   RANGE = 'V';
   UPLO = 'L';
   LDA  = n;
   ABSTOL = 0.0;  /* we use abstol = 0, since some lapack return an error otherwise */
   LDZ = n;
   M = -1;
   INFO = 0LL;

   /* interval of allowed values */
   VL = -1e30;
   VU = -tol;

   /* standard LAPACK workspace query, to get the amount of needed memory */
   LWORK = -1LL;
   LIWORK = -1LL;

   /* this computes the internally needed memory and returns this as (the first entries of [the 1x1 arrays]) WSIZE and WISIZE */
   F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
      &N, NULL, &LDA,
      &VL, &VU,
      NULL, NULL,
      &ABSTOL, &M, NULL, NULL,
      &LDZ, NULL, &WSIZE,
      &LWORK, &WISIZE, &LIWORK,
      &INFO);

   /* for some reason this code seems to be called with INFO=0 within UG */
   if ( convertToInt(INFO) != 0 )
   {
      SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d.\n", convertToInt(INFO));
      return SCIP_ERROR;
   }

   /* allocate workspace */
   LWORK = SCIP_RealTOINT(WSIZE);
   LIWORK = WISIZE;

   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WORK, (int) LWORK) );
   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &IWORK, (int) LIWORK) );
   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &ISUPPZ, (int) 2 * N) );

   /* call the function */
   F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
      &N, A, &LDA,
      &VL, &VU,
      NULL, NULL,
      &ABSTOL, &M, eigenvalues, eigenvectors,
      &LDZ, ISUPPZ, WORK,
      &LWORK, IWORK, &LIWORK,
      &INFO);

   /* free memory */
   BMSfreeBufferMemoryArray(bufmem, &ISUPPZ);
   BMSfreeBufferMemoryArray(bufmem, &IWORK);
   BMSfreeBufferMemoryArray(bufmem, &WORK);

   if ( convertToInt(INFO) != 0 )
   {
      SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d.\n", convertToInt(INFO));
      return SCIP_ERROR;
   }

   *neigenvalues = convertToInt(M);

   return SCIP_OKAY;
}


SCIP_RETCODE SCIPlapackComputeEigenvectorDecomposition(
   BMS_BUFMEM*           bufmem,             
   int                   n,                  
   SCIP_Real*            A,                  
   SCIP_Real*            eigenvalues,        
   SCIP_Real*            eigenvectors        
   )
{
   LAPACKINTTYPE* ISUPPZ;
   LAPACKINTTYPE* IWORK;
   LAPACKINTTYPE N;
   LAPACKINTTYPE INFO;
   LAPACKINTTYPE LDA;
   LAPACKINTTYPE LWORK;
   LAPACKINTTYPE LIWORK;
   LAPACKINTTYPE M;
   LAPACKINTTYPE LDZ;
   LAPACKINTTYPE WISIZE;
   SCIP_Real* WORK;
   SCIP_Real ABSTOL;
   SCIP_Real WSIZE;
   SCIP_Real VL;
   SCIP_Real VU;
   char JOBZ;
   char RANGE;
   char UPLO;

   assert( bufmem != NULL );
   assert( n > 0 );
   assert( n < INT_MAX );
   assert( A != NULL );
   assert( eigenvalues != NULL );
   assert( eigenvectors != NULL );

   N = n;
   JOBZ = 'V';
   RANGE = 'A';
   UPLO = 'L';
   LDA  = n;
   ABSTOL = 0.0;  /* we use abstol = 0, since some lapack return an error otherwise */
   M = n;
   LDZ = n;
   VL = -1e20;
   VU = 1e20;
   INFO = 0LL;

   /* standard LAPACK workspace query, to get the amount of needed memory */
   LWORK = -1LL;
   LIWORK = -1LL;

   /* this computes the internally needed memory and returns this as (the first entries of [the 1x1 arrays]) WSIZE and WISIZE */
   F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
      &N, NULL, &LDA,
      NULL, NULL,
      NULL, NULL,
      &ABSTOL, &M, NULL, NULL,
      &LDZ, NULL, &WSIZE,
      &LWORK, &WISIZE, &LIWORK,
      &INFO);

   if ( convertToInt(INFO) != 0 )
   {
      SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d.\n", convertToInt(INFO));
      return SCIP_ERROR;
   }

   /* allocate workspace */
   LWORK = SCIP_RealTOINT(WSIZE);
   LIWORK = WISIZE;

   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WORK, (int) LWORK) );
   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &IWORK, (int) LIWORK) );
   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &ISUPPZ, (int) 2 * N) );

   /* call the function */
   F77_FUNC(dsyevr, DSYEVR)( &JOBZ, &RANGE, &UPLO,
      &N, A, &LDA,
      &VL, &VU,
      NULL, NULL,
      &ABSTOL, &M, eigenvalues, eigenvectors,
      &LDZ, ISUPPZ, WORK,
      &LWORK, IWORK, &LIWORK,
      &INFO);

   /* free memory */
   BMSfreeBufferMemoryArray(bufmem, &ISUPPZ);
   BMSfreeBufferMemoryArray(bufmem, &IWORK);
   BMSfreeBufferMemoryArray(bufmem, &WORK);

   if ( convertToInt(INFO) != 0 )
   {
      SCIPerrorMessage("There was an error when calling DSYEVR. INFO = %d.\n", convertToInt(INFO));
      return SCIP_ERROR;
   }

   return SCIP_OKAY;
}


SCIP_RETCODE SCIPlapackMatrixVectorMult(
   int                   nrows,              
   int                   ncols,              
   SCIP_Real*            matrix,             
   SCIP_Real*            vector,             
   SCIP_Real*            result              
   )
{
   LAPACKINTTYPE M;
   LAPACKINTTYPE N;
   LAPACKINTTYPE LDA;
   LAPACKINTTYPE INCX;
   LAPACKINTTYPE INCY;
   SCIP_Real* A;
   SCIP_Real* X;
   SCIP_Real* Y;
   SCIP_Real BETA;
   SCIP_Real ALPHA;
   char TRANS;

   assert( nrows > 0 );
   assert( nrows < INT_MAX );
   assert( ncols > 0 );
   assert( ncols < INT_MAX );
   assert( matrix != NULL );
   assert( vector != NULL );
   assert( result != NULL );

   TRANS = 'N';
   M = nrows;
   N = ncols;
   ALPHA = 1.0;
   A = matrix;
   LDA = nrows;
   X = vector;
   INCX = 1;
   BETA = 0.0;
   Y = result;
   INCY = 1;

   F77_FUNC(dgemv, DGEMV)(&TRANS, &M, &N, &ALPHA, A, &LDA, X, &INCX, &BETA, Y, &INCY);

   return SCIP_OKAY;
}


SCIP_RETCODE SCIPlapackMatrixMatrixMult(
   int                   nrowsA,             
   int                   ncolsA,             
   SCIP_Real*            matrixA,            
   SCIP_Bool             transposeA,         
   int                   nrowsB,             
   int                   ncolsB,             
   SCIP_Real*            matrixB,            
   SCIP_Bool             transposeB,         
   SCIP_Real*            result              
   )
{
   LAPACKINTTYPE M;
   LAPACKINTTYPE N;
   LAPACKINTTYPE K;
   LAPACKINTTYPE LDA;
   LAPACKINTTYPE LDB;
   LAPACKINTTYPE LDC;
   SCIP_Real ALPHA;
   SCIP_Real BETA;
   char TRANSA;
   char TRANSB;

   assert( nrowsA > 0 );
   assert( nrowsA < INT_MAX );
   assert( ncolsA > 0 );
   assert( ncolsA < INT_MAX );
   assert( nrowsB > 0 );
   assert( nrowsB < INT_MAX );
   assert( ncolsB > 0 );
   assert( ncolsB < INT_MAX );
   assert( matrixA != NULL );
   assert( matrixB != NULL );
   assert( result != NULL );

   assert( (transposeA && transposeB && (nrowsA == ncolsB)) || (transposeA && !transposeB && (nrowsA == nrowsB))
      || (!transposeA && transposeB && (ncolsA == ncolsB)) || (!transposeA && !transposeB && (ncolsA == nrowsB)) );

   TRANSA = transposeA ? 'T' : 'N';
   TRANSB = transposeB ? 'T' : 'N';
   M = transposeA ? ncolsA : nrowsA;
   N = transposeB ? nrowsB : ncolsB;
   K = transposeA ? nrowsA : ncolsA;
   ALPHA = 1.0;
   LDA = transposeA ? K : M;
   LDB = transposeB ? N : K;
   BETA = 0.0;
   LDC = M;

   F77_FUNC(dgemm, DGEMM)(&TRANSA, &TRANSB, &M, &N, &K, &ALPHA, matrixA, &LDA, matrixB, &LDB, &BETA, result, &LDC);

   return SCIP_OKAY;
}


SCIP_RETCODE SCIPlapackLinearSolve(
   BMS_BUFMEM*           bufmem,             
   int                   m,                  
   int                   n,                  
   SCIP_Real*            A,                  
   SCIP_Real*            b,                  
   SCIP_Real*            x                   
   )
{
   LAPACKINTTYPE* IWORK;
   LAPACKINTTYPE M;
   LAPACKINTTYPE N;
   LAPACKINTTYPE NRHS;
   LAPACKINTTYPE LDA;
   LAPACKINTTYPE LDB;
   LAPACKINTTYPE RANK;
   LAPACKINTTYPE LWORK;
   LAPACKINTTYPE LIWORK;
   LAPACKINTTYPE INFO;
   LAPACKINTTYPE WISIZE;
   SCIP_Real* WORK;
   SCIP_Real* S;
   SCIP_Real RCOND;
   SCIP_Real WSIZE;
   int i;

   assert( bufmem != NULL );
   assert( m > 0 );
   assert( m < INT_MAX );
   assert( n > 0 );
   assert( n < INT_MAX );
   assert( A != NULL );
   assert( b != NULL );
   assert( x != NULL );

   M = m;
   N = n;
   NRHS = 1;
   LDA = m;
   LDB = MAX(M,N);
   RCOND = 0.0;
   INFO = 0LL;

   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &S, MIN(m,n)) );

   /* standard LAPACK workspace query, to get the amount of needed memory */
   LWORK = -1LL;

   /* this computes the internally needed memory and returns this as (the first entry of [the 1x1 array]) WSIZE */
   F77_FUNC(dgelsd, DGELSD)( &M, &N, &NRHS,
      NULL, &LDA, NULL, &LDB, NULL,
      &RCOND, &RANK, &WSIZE, &LWORK,
      &WISIZE, &INFO);

   if ( convertToInt(INFO) != 0 )
   {
      SCIPerrorMessage("There was an error when calling DGELSD. INFO = %d\n", convertToInt(INFO));
      return SCIP_ERROR;
   }

   /* allocate workspace */
   LWORK = SCIP_RealTOINT(WSIZE);
   LIWORK = WISIZE;

   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &WORK, (int) LWORK) );
   BMS_CALL( BMSallocBufferMemoryArray(bufmem, &IWORK, (int) LIWORK) );

   /* call the function */
   F77_FUNC(dgelsd, DGELSD)( &M, &N, &NRHS,
      A, &LDA, b, &LDB, S, &RCOND, &RANK,
      WORK, &LWORK, IWORK, &INFO);

#ifdef PRINTMATRICES
   {
      SCIP_Real residual = 0.0;

      printf("LWORK = %d\n", LWORK);
      printf("LIWORK = %d\n", LIWORK);
      printf("A has size (%d,%d), is of rank %d\n", M, N, RANK);
      printf("Minimum l2-norm solution of linear equation system:\n");

      for (i = 0; i < n; ++i)
         printf("(%d, %f)   ", i, b[i]);

      for (i = n; i < m; ++i)
         residual += b[i] * b[i];
      printf("\n");

      printf("Residual sum-of-squares for the solution is %f\n", residual);
   }
#endif

   /* free memory */
   BMSfreeBufferMemoryArray(bufmem, &IWORK);
   BMSfreeBufferMemoryArray(bufmem, &WORK);
   BMSfreeBufferMemoryArray(bufmem, &S);

   if ( convertToInt(INFO) != 0 )
   {
      SCIPerrorMessage("There was an error when calling DGELSD. INFO = %d\n", convertToInt(INFO));
      return SCIP_ERROR;
   }

   /* LAPACK overwrites the right-hand side with the result */
   for (i = 0; i < n; ++i)
      x[i] = b[i];

   return SCIP_OKAY;
}