[Insight-users] MultiplyMatrixVector in LinearSystemWrapperVNL

[Jeffrey Duda]

Hi,
Its hard to say exactly why the matrix vector multiplication isn't working without seeing the rest of your code.  I've included a quick example of using the LinearSystemWrapperVNL to do the operation that you appear to be attempting.  The output of this example is shown below the code.  Hopefully you can compare your code to this and figure out where things are going awry.

-jeffrey duda

---------------------------------------

#include "itkFEMLinearSystemWrapperVNL.h"
#include <iostream>

using namespace itk::fem;
using namespace std;

/* Testing for linear system wrappers */
int main( void ) 
{

  /* loop vars for printing */
  unsigned int i;
  unsigned int j;

  /* declare wrapper */
  LinearSystemWrapperVNL it;

  /* system parameters */
  unsigned int N = 5;
  unsigned int nMatrices =  1;
  unsigned int nVectors =   2;

  /* Set up the system */
  it.SetSystemOrder(N);
  it.SetNumberOfMatrices(nMatrices);
  it.SetNumberOfVectors(nVectors);  

  /* Initialize memory */
  for (i=0; i<nMatrices; i++) 
  {
    it.InitializeMatrix(i);
  }
  for (i=0; i<nVectors; i++)
  {
    it.InitializeVector(i);
  }

  /*     matrix 0
   * |11  0  0 14 15|
   * | 0 22  0  0  0|
   * | 0  0 33  0  0|
   * |14  0  0 44 45|
   * |15  0  0 45 55|
   */
  it.SetMatrixValue(0,0,11,0); it.SetMatrixValue(0,3,14,0); it.SetMatrixValue(0,4,15,0);
  it.SetMatrixValue(1,1,22,0); 
  it.SetMatrixValue(2,2,33,0); 
  it.SetMatrixValue(3,0,14,0); it.SetMatrixValue(3,3,44,0); it.SetMatrixValue(3,4,45,0);
  it.SetMatrixValue(4,0,15,0); it.SetMatrixValue(4,3,45,0); it.SetMatrixValue(4,4,55,0);

  /* print matrix 0 */
  cout << "Matrix 0" << endl;
  for(i=0; i<N; i++) 
  {
    for (j=0; j<N; j++) 
    {
      cout << it.GetMatrixValue(i,j,0) << " ";
    }
    cout << endl;
  }
  cout << endl;

  /* Vector 0 = [1 2 3 4 5] */
  it.SetVectorValue(0,1,0);
  it.SetVectorValue(1,2,0);
  it.SetVectorValue(2,3,0);
  it.SetVectorValue(3,4,0);
  it.SetVectorValue(4,5,0);

  /* print Vector 0 */
  cout << "Vector 0" << endl;
  for (i=0; i<N; i++) 
  {
    cout << it.GetVectorValue(i,0) << " ";
  }
  cout << endl << endl;

  /* vector 1 = matrix 0 * vector 0 */
  cout << "Vector 1 =  Matrix 0 * Vector 0" << endl;
  it.MultiplyMatrixVector(1, 0, 0);

  /* print Vector 1 */
  for (i=0; i<N; i++) 
  {
    cout << it.GetVectorValue(i,1) << " ";
  }
  cout << endl << endl;

  /* destroy matrix,vector,solution */
  it.DestroyMatrix(0);
  it.DestroyVector(0);
  it.DestroyVector(1);

  return 0;
 
}

-------- output of above code -----------

Matrix 0
11 0 0 14 15 
0 22 0 0 0 
0 0 33 0 0 
14 0 0 44 45 
15 0 0 45 55 

Vector 0
1 2 3 4 5 

Vector 1 =  Matrix 0 * Vector 0
142 44 99 415 470 

-----Original Message-----
> From: insight-users-admin@public.kitware.com 
> [mailto:insight-users-admin@public.kitware.com] On Behalf Of 
> du Bois d'Aische
> Sent: Friday, April 04, 2003 1:40 PM
> To: insight-users@public.kitware.com
> Subject: [Insight-users] MultiplyMatrixVector in 
> LinearSystemWrapperVNL
> 
> 
> 
>   Hi,
> 
>       I try to use to method MultiplyMatrixVector in 
> LinearSystemWrapperVNL. This method requires an unsigned int 
> resultVectorIndex, an unsigned int matrixIndex and an 
> unsigned int vectorIndex I have initialized the vectors and 
> the matrix 
>     m_ls->SetSystemOrder(length);
>     m_ls->SetNumberOfMatrices(1);
>     m_ls->InitializeMatrix(Ma);
>     m_ls->SetNumberOfVectors(2);
>     m_ls->InitializeVector(Ve);
>     m_ls->InitializeVector(Veresu);
>  where m_ls is the Pointer to the LinearSystemWrapperVNL; Ma, 
> Ve and Veresu are unsigned int. Then I fill the matrix Ma and 
> the vector Ve. I display them and they look ok, 
>    but the result of m_ls->MultiplyMatrixVector(Veresu,Ma,Ve) 
> is always a vector of the size 'length' full of 0;
> 
>            Has someone an idea why the result could be always 0 ?
> 
> 
>                   thanks,
> 
>                        aloys 
> _______________________________________________
> Insight-users mailing list
> Insight-users@public.kitware.com 
> http://public.kitware.com/mailman/listinfo/ins> ight-users