Added function to perform QR matrix decomposition

ayendiki · Apr 21, 2020 · 3f8ab07 · 3f8ab07
1 parent be5233c
commit 3f8ab07
Show file tree

Hide file tree

Showing 4 changed files with 55 additions and 0 deletions.
diff --git a/include/matrix.h b/include/matrix.h
@@ -85,6 +85,7 @@ int     MatrixIsZero(MATRIX *m) ;
 int     MatrixIsIdentity(MATRIX *m) ;
 MATRIX  *MatrixReshape(MATRIX *m_src, MATRIX *m_dst, int rows, int cols) ;
 int     MatrixCheck(MATRIX *m) ;
+int     MatrixQRdecomposition(const MATRIX *iMatrix, MATRIX *oQ, MATRIX *oR);
 MATRIX  *MatrixInverse( const MATRIX *mIn, MATRIX *mOut) ;
 MATRIX  *MatrixPseudoInverse(MATRIX *m, MATRIX *m_pseudo_inv) ;
 MATRIX  *MatrixSVDPseudoInverse(MATRIX *m, MATRIX *m_pseudo_inv) ;

diff --git a/include/numerics.h b/include/numerics.h
@@ -38,6 +38,8 @@
   // this is only called from the matrix class
   int OpenSvdcmp( MATRIX *a, VECTOR *w, MATRIX *v );
 
+  int OpenQRdecomposition(const MATRIX *iMatrix, MATRIX *oQ, MATRIX *oR);
+
   float OpenMatrixDeterminant( MATRIX *matrix );
 
   int OpenLUMatrixInverse( MATRIX *matrix, MATRIX *inverse );

diff --git a/utils/matrix.cpp b/utils/matrix.cpp
@@ -175,6 +175,12 @@ MATRIX *MatrixCopy(const MATRIX *mIn, MATRIX *mOut)
   return (mOut);
 }
 
+int MatrixQRdecomposition(const MATRIX *iMatrix, MATRIX *oQ, MATRIX *oR)
+{
+  int r = OpenQRdecomposition(iMatrix,oQ,oR);
+  return(r);
+}
+
 MATRIX *MatrixInverse(const MATRIX *mIn, MATRIX *mOut)
 {
   // float **a, **y;
@@ -4707,3 +4713,4 @@ MATRIX *MatrixGlmFit(MATRIX *y, MATRIX *X, double *pRVar, MATRIX *beta)
   MatrixFree(&res);
   return(beta);
 }
+
diff --git a/utils/numerics.cpp b/utils/numerics.cpp
@@ -1145,6 +1145,51 @@ int OpenLUMatrixInverse(MATRIX *iMatrix, MATRIX *oInverse)
   return errorCode;
 }
 
+/*!
+  \fn int OpenQRdecomposition(const MATRIX *iMatrix, MATRIX *oQ, MATRIX *oR)
+  \brief Performs QR decomposition on the input matrix. The output
+  matricies must be allocated already (Q is nrowsxnows where nrows is
+  the number of rows in the input matrix, and R is same size as the
+  input matrix). 
+ */
+int OpenQRdecomposition(const MATRIX *iMatrix, MATRIX *oQ, MATRIX *oR)
+{
+  if(iMatrix==NULL){
+    printf("QRdecomposition(): input matrix is NULL\n");
+    return(1);
+  }
+  if(oQ==NULL){
+    printf("QRdecomposition(): output Q matrix is NULL\n");
+    return(1);
+  }
+  if(oR==NULL){
+    printf("QRdecomposition(): output R matrix is NULL\n");
+    return(1);
+  }
+  if(oQ->rows != iMatrix->rows){
+    printf("QRdecomposition(): Q rows != input Matrix rows\n");
+    return(1);
+  }
+  if(oQ->cols != oQ->rows){
+    printf("QRdecomposition(): Q cols != Q rows\n");
+    return(1);
+  }
+  if(oR->rows != iMatrix->rows){
+    printf("QRdecomposition(): R rows != input Matrix rows\n");
+    return(1);
+  }
+  if(oR->cols != iMatrix->cols){
+    printf("QRdecomposition(): R cols != input Matrix cols\n");
+    return(1);
+  }
+
+  vnl_matrix< float >  vnlMatrix(iMatrix->data, iMatrix->rows, iMatrix->cols);
+  vnl_qr<float> qr( vnlMatrix );
+  qr.Q().copy_out(oQ->data);
+  qr.R().copy_out(oR->data);
+  return(0);
+}
+
 /**
  * Returns the determinant of matrix.
  * @param iMatrix