triangularize a matrix in one go with drotg & dlasr

[bravegag]

Hello,

I have been trying to use Givens rotations via a combination of drotg and dlasr but have not been able to get it right. Can anyone please help me out. Please find below the self contained code of what I am trying to acomplish.

Many thanks in advance,
Best regards,
Giovanni

Code: Select all

/*
 * test.cc
 *
 *  Created on: May 31, 2012
 *      Author: Giovanni Azua
 */
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#include <mkl_lapack.h>
#include <mkl_lapacke.h>
#include <mkl_cblas.h>

static const int ROWS = 3;
static const int COLS = 3;
static double matrix[ROWS][COLS];

static void print_matrix() {
   int rows = ROWS;
   int cols = COLS;
   printf("matrix of dimensions %dx%d: \n", rows, cols);
   for (int i=0; i < rows; ++i) {
      for (int j=0; j < cols; ++j) {
         printf("%f ", matrix[j][i]);
      }
      printf("\n");
   }
}

int main(int argc, char** argv) {
   // create a simple matrix in "column major" with 1's in the
   // upper triangule and 2s in the band below the diagonal
   // now I want to get rid of the 2s using Givens rotations
   matrix[0][0] = 1; matrix[1][0] = 1; matrix[2][0] = 1;
   matrix[0][1] = 2; matrix[1][1] = 1; matrix[2][1] = 1;
   matrix[0][2] = 0; matrix[1][2] = 2; matrix[2][2] = 1;
   print_matrix();

   double c[COLS];
   double s[COLS];

   // create the rotations
   for (int j = 0; j < (COLS - 1); j++) {
      double a = matrix[j][j];
      double b = matrix[j + 1][j];
      cblas_drotg(&a, &b, &c[j], &s[j]);
   }

   // apply the rotations along the diagonal
   char side      = 'L';
   char pivot     = 'V';
   char direct    = 'F';
   lapack_int m   = ROWS;
   lapack_int n   = COLS;
   lapack_int lda = ROWS;
   dlasr(&side, &pivot, &direct, &m, &n, c, s, (double*) matrix, &lda);
   print_matrix();

   return 0;
}


[Julien Langou]

Hello,

I glanced at the code quickly. I assume that you take a Hessenberg matrix and then you want to triangularize it through Givens rotations.
If this is the case, the rotation generation and the rotation application needs to be interlaced.
You need to create the first rotation, then apply it, then create the second, then apply it, etc.
Because the first rotation affects the second and the third row and so you need to apply it before creating the second rotation which needs elements in the third row.
Not sure this helps though. As I said as just glanced at the code.

Julien.

[bravegag]

Hi Julien,

Thank you. Indeed, it is a upper Hessenberg matrix. I was actually aiming to do the triangulation in a "BLAS level 2" fashion that's why I was trying to use dlasr. Do you have example code for using this function dlasr?

From the top of your head can you point me to one example of this use-case of trying to triangularize in "one go" and using Given rotations a matrix in upper Hessenberg form i.e. when one column in the middle has been deleted?

Many TIA,
Best regards,
Giovanni

[bravegag]

Never mind, figured it out already. Please see below the working example self-contained C++ code:

Code: Select all

/*
 * test.cc
 *
 *  Created on: May 31, 2012
 *      Author: Giovanni Azua
 */
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#include <mkl_lapack.h>
#include <mkl_lapacke.h>
#include <mkl_cblas.h>

static const int ROWS = 3;
static const int COLS = 3;
static double matrix[ROWS][COLS];

static void print_matrix() {
   int rows = ROWS;
   int cols = COLS;
   printf("matrix of dimensions %dx%d: \n", rows, cols);
   for (int i=0; i < rows; ++i) {
      for (int j=0; j < cols; ++j) {
         printf("%f ", matrix[j][i]);
      }
      printf("\n");
   }
}

int main(int argc, char** argv) {

   // create a simple matrix in "column major" with 1's in the
   // upper triangule and 2s in the band below the diagonal
   // now I want to get rid of the 2s using Givens rotations

   matrix[0][0] = 1; matrix[1][0] = 1; matrix[2][0] = 1;
   matrix[0][1] = 2; matrix[1][1] = 1; matrix[2][1] = 1;
   matrix[0][2] = 0; matrix[1][2] = 2; matrix[2][2] = 1;
   print_matrix();

   // Givens's sins and cos
   double cc[COLS];
   double ss[COLS];

   // create the rotations
   for (int j = 0; j < COLS; j++) {
      double a = matrix[j][j];
      double b = matrix[j][j + 1];
      double c = 0.0;
      double s = 0.0;
      double r = 0.0;
      dlartg(&a, &b, &c, &s, &r);
      // could also use cblas_drotg
      //cblas_drotg(&a, &b, &c, &s);

      for (int jj = j; jj < COLS; ++jj) {
         cc[jj] = c;
         ss[jj] = s;
      }

      // apply the rotations to trailing columns on the right
      char side      = 'L';
      char pivot     = 'V';
      char direct    = 'F';
      lapack_int m   = ROWS - j;
      lapack_int n   = COLS - j;
      lapack_int lda = ROWS;
      dlasr(&side, &pivot, &direct, &m, &n, &cc[j], &ss[j], ((double*) matrix) + j*ROWS + j, &lda);

      print_matrix();
   }

   return 0;
}