/*  ldumat.c    CCMATH mathematics library source code.
 *
 *  Copyright (C)  2000   Daniel A. Atkinson    All rights reserved.
 *  This code may be redistributed under the terms of the GNU library
 *  public license (LGPL). ( See the lgpl.license file for details.)
 * ------------------------------------------------------------------------
 */
#include <stdlib.h>
void ldumat(double *a,double *u,int m,int n)
{ double *p0,*q0,*p,*q,*w;
  int i,j,k,mm;
  double s,h;
  w=(double *)calloc(m,sizeof(double));
  for(i=0,mm=m*m,q=u; i<mm ;++i) *q++ =0.;
  p0=a+n*n-1; q0=u+m*m-1; mm=m-n; i=n-1;
  for(j=0; j<mm ;++j,q0-=m+1) *q0=1.;
  if(mm==0){ p0-=n+1; *q0=1.; q0-=m+1; --i; ++mm;}
  for(; i>=0 ;--i,++mm,p0-=n+1,q0-=m+1){
    if(*p0!=0.){
      for(j=0,p=p0+n,h=1.; j<mm ;p+=n) w[j++]= *p;
      h= *p0; *q0=1.-h;
      for(j=0,q=q0+m; j<mm ;q+=m) *q= -h*w[j++];
      for(k=i+1,q=q0+1; k<m ;++k){
	for(j=0,p=q+m,s=0.; j<mm ;p+=m) s+=w[j++]* *p;
	s*=h;
	for(j=0,p=q+m; j<mm ;p+=m) *p-=s*w[j++];
        *q++ = -s;
       }
     }
    else{
      *q0=1.;
      for(j=0,p=q0+1,q=q0+m; j<mm ;++j,q+=m) *q= *p++ =0.;
     }
   }
  free(w);
}