← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/control/wrref.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/control/wrref.f
 
1
C/MEMBR ADD NAME=WRREF,SSI=0
 
2
c     Copyright INRIA
 
3
      subroutine wrref(ar,ai,lda,m,n,eps)
 
4
c!but
 
5
c     wrref calcule la forme echelon d'une matrice a coeff complexes
 
6
c!liste d'appel
 
7
c
 
8
c     subroutine wrref(ar,ai,lda,m,n,eps)
 
9
c     double precision ar(lda,n),ai(lda,n),eps
 
10
c     integer lda,m,n
 
11
c
 
12
c     ar,ai : tableaux contenant a l'appel les parties reelles et
 
13
c        complexes  de la matrice dont on veut determiner la forme
 
14
c        echelon, apres execution a contient la forme echelon
 
15
c     lda : nombre de ligne du tableau a dans le programme appelant
 
16
c     m : nombre de ligne de la matrice a
 
17
c     n : nombre de colonnes de a matrice a
 
18
c     eps : tolerance. les reels  inferieurs a 2*max(m,n)*eps*n1(a),
 
19
c           ou n1(a) est la norme 1 de a ,sont consideres comme nuls
 
20
c
 
21
c     si l'on veut la transformation appliquee,appeler wrref avec
 
22
c     la matrice obtenue en concatenant l'identite a la matrice a
 
23
c     en rajoutant des colonnes.
 
24
c!sous programmes appeles
 
25
c     iwamax wcopy wswap wscal wasum waxpy (blas.extensions)
 
26
c     dset (blas)
 
27
c     dble abs max (fortran)
 
28
c!
 
29
c!
 
30
      double precision ar(lda,*),ai(lda,*),eps,tol,tr,ti,wasum
 
31
      tol = 0.0d+0
 
32
      do 10 j = 1, n
 
33
         tol = max(tol,wasum(m,ar(1,j),ai(1,j),1))
 
34
   10 continue
 
35
      tol = eps*dble(2*max(m,n))*tol
 
36
      k = 1
 
37
      l = 1
 
38
   20 if (k.gt.m .or. l.gt.n) return
 
39
      i = iwamax(m-k+1,ar(k,l),ai(k,l),1) + k-1
 
40
      if (abs(ar(i,l))+abs(ai(i,l)) .gt. tol) go to 30
 
41
      call dset(m-k+1,0.0d+0,ar(k,l),1)
 
42
      call dset(m-k+1,0.0d+0,ai(k,l),1)
 
43
         l = l+1
 
44
         go to 20
 
45
   30 call wswap(n-l+1,ar(i,l),ai(i,l),lda,ar(k,l),ai(k,l),lda)
 
46
      call wdiv(1.0d+0,0.0d+0,ar(k,l),ai(k,l),tr,ti)
 
47
      call wscal(n-l+1,tr,ti,ar(k,l),ai(k,l),lda)
 
48
      ar(k,l) = 1.0d+0
 
49
      ai(k,l) = 0.0d+0
 
50
      do 40 i = 1, m
 
51
         tr = -ar(i,l)
 
52
         ti = -ai(i,l)
 
53
         if (i .ne. k) call waxpy(n-l+1,tr,ti,
 
54
     $                 ar(k,l),ai(k,l),lda,ar(i,l),ai(i,l),lda)
 
55
   40 continue
 
56
      k = k+1
 
57
      l = l+1
 
58
      go to 20
 
59
      end