void umiset(int n, double *a, int lda)
#define A(I,J) a[(I) + (J)* ( lda)]
{
  int i, j;

  for (i=0; i<=n-1; i++) {
     for (j=0; j<=n-1; j++) 
        A(i,j)=0;
      A(i,i)=1;
  } 
}   

void umoset(int m, int n, double *a, int lda)
#define A(I,J) a[(I) + (J)* ( lda)]
{
  int i, j;
  for (i=0; i<=m-1; i++) 
     for (j=0; j<=n-1; j++) 
        A(i,j)=0;
}