source: branches/GRISLIv3/SOURCES/New-remplimat/eq_ellipt_sgbsv_mod-0.2.f90 @ 406

Last change on this file since 406 was 401, checked in by aquiquet, 16 months ago

Cleaning branch: use only statements in New-remplimat
File size: 21.1 KB
Line 
1module eq_ellip_sgbsv_mod
2
3!   module pour la resolution des equations elliptiques
4!   de l'ice shelf et ice stream en utilisant une resolution
5!   de systeme bande  sgbsv (lapack)
6
7!   voir le rappel des argument de SGBSV en fin de fichier
8
9! declaration des variables
10!$ USE OMP_LIB
11use module3D_phy,    only: nlmin
12use remplimat_declar,only: nflmax
13
14implicit none
15
16integer,parameter :: klmax=2*nlmin           ! largeur max demi-bande inferieure
17integer,parameter :: kumax=2*nlmin           ! largeur max demi-bande superieure
18integer,parameter :: ldbmx=2*klmax+kumax+1   ! largeur de bande maximum (taille de la matrice) >= 2kl+ku+1
19
20integer :: lbande        ! pour appel sgbsv (largeur de bande) = kl+ku+1 (ne pas confondre avec ldbmx)
21integer :: ldb           ! pour appel sgbsv (largeur de bande+ travail) = 2*kl+ku+1
22integer :: ldb_old       ! pour garder la valeur du tour precedent (allocation de mmat)
23integer :: nrhs = 1      ! en général 1
24
25
26integer :: ku            ! largeur demi-bande superieure effective
27integer :: kl            ! largeur demi-bande inferieure effective
28integer :: nblig         ! nombre d'equations traitees effectivement
29integer :: nblig_old     ! pour garder la valeur du tour precedent (allocation de mmat)
30
31integer :: ldiag         ! compteur du numero de ligne
32 
33integer :: iii,jjj       ! position du noeud de l'echec de sgbsv
34integer :: ifail_sgbsv   ! recuperation d'erreur
35
36
37
38! tableaux pour passer a sgbsv
39
40real,dimension(nflmax,1)         :: bdr       ! vecteur membre de droite de l'equation
41integer,dimension(nflmax)        :: ipiv      ! pivot
42
43
44! seul mmat est allocatable parce que la premiere dimension est variable
45
46real,dimension(:,:), allocatable ::  mmat     ! AB, matrice appelee par sgbsv contenant les diagonales
47                                              ! sous forme "couchees" sera dimensionnee (ldab,nblig)
48
49! tableaux de travail
50integer,dimension(nflmax) :: lu_band          ! largeur de bande
51integer,dimension(nflmax) :: ll_band          ! largeur de bande inferieure
52
53
54
55contains
56
57
58!----------------------------
59subroutine  redim_matrice
60
61! allocation dynamique de mmat
62
63implicit none
64
65integer ::  err
66
67if (allocated(mmat)) then
68   deallocate(mmat,stat=err)
69     if (err/=0) then
70        print *,"Erreur à la desallocation de mmat",err
71        stop
72     end if
73end if
74
75
76allocate(mmat(ldb,nblig),stat=err)
77if (err/=0) then
78   print *,"erreur a l'allocation du tableau mmat",err
79   print *,"ldb,nblig",ldb,nblig
80   stop
81end if
82
83
84!write(6,*) "allocation mmat ok, ldb,nblig",ldb,nblig
85
86end subroutine redim_matrice
87
88! -------------------------------------------------------------------------------
89
90subroutine resol_ellipt(nx1,nx2,ny1,ny2,uxprec,uyprec,uxnew,uynew,ifail)
91
92!  prepare la matrice mmat (version "couchee" de L2)
93!  resoud l'equation elliptique par un appel a sgbsv
94!  renvoie les nouvelles vitesses
95!
96!                                                                   
97!      uxprex(n1,n2)                                                 
98!      uyprec(n1,n2) vitesses de l'iteration precedente             
99!                                                                   
100!      uxnew(n1,n2)                                                   
101!      uynew(n1,n2) uynew resultat de cette iteration
102!
103!---------------------------------------------------------------------------
104
105use module3D_phy, only: debug_3D
106use remplimat_declar, only: ok_umat,ok_vmat,ligu_l2,ligv_l2
107
108implicit none
109
110integer,intent(in) :: nx1             ! bornes du domaine en x (noeuds majeurs)
111integer,intent(in) :: nx2
112integer,intent(in) :: ny1             ! bornes du domaine en y (noeuds majoeurs)
113integer,intent(in) :: ny2
114
115integer :: ifail                      ! pour les rapports d'erreur
116
117real, dimension(nx1:nx2,ny1:ny2),intent(in)  :: uxprec    ! vitesse en entree routine
118real, dimension(nx1:nx2,ny1:ny2),intent(in)  :: uyprec    ! vitesse en entree routine
119
120real, dimension(nx1:nx2,ny1:ny2),intent(out) :: uxnew     ! vitesse en sortie de la routine
121real, dimension(nx1:nx2,ny1:ny2),intent(out) :: uynew     ! vitesse en sortie de la routine
122
123integer i,j ! loop integers
124
125nrhs = 1
126nblig_old=0
127ldb_old=0
128
129call prepare_mat(nx1,nx2,ny1,ny2)         ! prepare la matrice pour appel a SGBSV
130
131! quelques outils de debug ci-dessous
132!---------------------------------------
133
134! write(6,'(a20,4(i0,2x))') 'kl, ku ,ldb, nblig   ',kl, ku ,ldb, nblig
135!---------------------------------------------------------
136
137! routine qui permet de visualiser les mmatrices L2 et mmat
138! call graphique_L2(kl,ku,nx1,nx2,ny1,ny2,imx(nx1:nx2,ny1:ny2),imy(nx1:nx2,ny1:ny2))   
139!---------------------------------------------------------
140
141! routine qui permet de tester l'appel sgbsv et de visualiser la matrice
142! call graph_sgbsv(nblig,kl,ku,nrhs,mmat,ldb,ipiv,bdr,nblig,ifail_sgbsv)
143!---------------------------------------------------------
144
145
146
147! appel a la resolution du systeme (voir doc a la fin du fichier)
148
149
150if (nblig.ne.0) then
151   call sgbsv(nblig,kl,ku,nrhs,mmat,ldb,ipiv,bdr,nblig,ifail_sgbsv)
152
153   if (ifail_sgbsv.eq.0) then                ! la resolution s'est bien passee
154      !   write(6,*)'-----------------------'
155      !   write(6,*) 'sortie ok de sgbsv'
156      !   write(6,*)'-----------------------'
157      !$OMP PARALLEL
158      !$OMP DO
159      do j = ny1,ny2
160         do i = nx1,nx2
161
162            if (ok_umat(i,j)) then
163               uxnew(i,j) = bdr(ligu_L2(i,j),1)
164            else
165               uxnew(i,j) = uxprec(i,j)
166            end if
167
168            if (ok_vmat(i,j)) then
169               uynew(i,j) = bdr(ligv_L2(i,j),1)
170            else
171               uynew(i,j) = uyprec(i,j)
172            end if
173         end do
174      end do
175      !$OMP END DO
176      !$OMP END PARALLEL
177
178   else
179
180      call erreur_sgbsv(ifail_sgbsv)
181      ifail=ifail_sgbsv
182   end if
183
184else                 ! nblig=0, pas de passage par sgbsv
185   !$OMP PARALLEL
186   !$OMP WORKSHARE
187   uxnew(nx1:nx2,ny1:ny2)=uxprec(nx1:nx2,ny1:ny2)
188   uynew(nx1:nx2,ny1:ny2)=uyprec(nx1:nx2,ny1:ny2)
189   !$OMP END WORKSHARE
190   !$OMP END PARALLEL
191end if
192
193debug_3D(nx1:nx2,ny1:ny2,37)=uxnew(nx1:nx2,ny1:ny2)
194debug_3D(nx1:nx2,ny1:ny2,38)=uynew(nx1:nx2,ny1:ny2)
195
196return
197end subroutine resol_ellipt
198
199!--------------------------------------------------------------------------------
200
201subroutine prepare_mat(nx1,nx2,ny1,ny2)         ! prepare la matrice pour appel a SGBSV
202
203use remplimat_declar, only: count_line,ligu_l2,ligv_l2,pos_ligu,pos_ligv,ok_umat,ok_vmat, &
204                            il,jl,ilmin,ilmax,jlmin,jlmax,tu,tv,su,sv,ghost_x,ghost_y,    &
205                            opposx,opposy,eps_col
206
207implicit none
208
209integer :: nx1,nx2        ! bornes du domaine traite
210integer :: ny1,ny2        ! bornes du domaine traite
211
212
213integer :: ligne          ! dans L2
214integer :: colonne        ! dans L2
215integer :: ll             ! dans le calcul de la largeur de bande
216
217! tableau de travail pour calculer la largeur de bande
218integer,dimension(20)     :: pos
219
220integer i,j ! loop integers
221
222! numerotation des lignes
223!--------------------------
224!$OMP PARALLEL
225!$OMP WORKSHARE
226ligu_L2(:,:) = -9999
227ligv_L2(:,:) = -9999
228pos_ligu(:,:) = -9999
229pos_ligv(:,:) = -9999
230!$OMP END WORKSHARE
231!$OMP END PARALLEL
232! Boucle sur les noeuds avec la dimension la plus grande a l'interieur
233
234count_line = 1                                 ! pour numeroter les lignes
235
236count: do j = ny1,ny2
237   do i = nx1,nx2
238
239      if (ok_umat(i,j)) then
240         ligu_L2(i,j) = count_line             ! ligne de U(i,j) dans L2
241
242         pos_ligu(count_line,1) = i            ! i,j d'un noeud U de ligne count_line
243         pos_ligu(count_line,2) = j
244
245         count_line = count_line+1
246
247
248                                       
249      end if
250 
251      if (ok_vmat(i,j)) then
252         ligv_L2(i,j) = count_line             ! ligne de V(i,j) dans L2
253
254         pos_ligv(count_line,1) = i            ! i,j d'un noeud V de ligne count_line
255         pos_ligv(count_line,2) = j
256
257         count_line = count_line+1
258
259      end if
260
261
262   end do
263end do count
264
265
266nblig = count_line-1
267
268
269! calcul de la largeur de bande. depend de la numerotation des lignes
270!--------------------------------------------------------------------
271!$OMP PARALLEL
272!$OMP WORKSHARE
273lu_band(:) = 0
274ll_band(:) = 0
275!$OMP END WORKSHARE
276!$OMP END PARALLEL
277
278larg_band: do j = ny1,ny2
279   do i = nx1,nx2
280
281      lband_U: if (ok_umat(i,j)) then         ! pour chaque ligne uij, on rempli pos avec les numeros de colonne
282
283!     ilmin et jlmin boucle U
284      ilmin = max(-1,nx1+1-i)    ! pour avoir i+il entre nx1+1 et nx2
285      ilmax = min(1,nx2-i)
286
287      jlmin = max(-1,ny1-j)      ! pour avoir j+jl entre ny1 et ny2
288      jlmax = min(1,ny2-j) 
289
290   
291         ll = 1 
292         ldiag = ligu_L2(i,j)            ! numero colonne de la diagonale 
293         pos(:) = ldiag             
294
295         do jl = jlmin , jlmax      ! balaye tous les coefficients de la colonne Uij de -1 a 1
296            do il = ilmin , ilmax   !                                                de -1 a 1
297
298               if (Tu(i,j,il,jl).ne.0.) then
299                  pos(ll) = ligu_L2(i+il,j+jl)
300                  ll = ll+1
301               end if
302               
303               if (Tv(i,j,il,jl).ne.0.) then
304                  pos(ll) = ligv_L2(i+il,j+jl)
305                  ll = ll+1
306               end if
307            end do
308         end do
309
310
311         pos(:) = pos(:)-ldiag    ! on soustrait le numero de colonne de celui de la diagonale
312         ll_band(ldiag) = -minval(pos)
313         lu_band(ldiag) = maxval(pos)
314
315      end if lband_U
316
317      lband_V: if (ok_vmat(i,j)) then         ! pour chaque ligne uij, on rempli pos avec les numeros de colonne
318
319!     ilmin et jlmin boucle V
320      ilmin = max(-1,nx1-i)    ! pour avoir i+il entre nx1 et nx2
321      ilmax = min(1,nx2-i)
322
323      jlmin = max(-1,ny1+1-j)    ! pour avoir j+jl entre ny1+1 et ny2
324      jlmax = min(1,ny2-j)     
325
326       
327         ll = 1 
328         ldiag = ligv_L2(i,j)            ! numero colonne de la diagonale 
329         pos(:) = ldiag             
330
331         do jl = jlmin , jlmax      ! balaye tous les coefficients de la colonne Uij de -1 a 1
332            do il = ilmin , ilmax   !                                                de -1 a 1
333
334               if (Su(i,j,il,jl).ne.0.) then
335                  pos(ll) = ligu_L2(i+il,j+jl)
336                  ll = ll+1
337               end if
338               
339               if (Sv(i,j,il,jl).ne.0) then
340                  pos(ll) = ligv_L2(i+il,j+jl)
341                  ll = ll+1
342               end if
343            end do
344         end do
345
346
347         pos(:) = pos(:)-ldiag    ! on soustrait le numero colonne de la diagonale
348         ll_band(ldiag) = -minval(pos)
349         lu_band(ldiag) = maxval(pos)
350
351      end if lband_V
352   end do
353end do larg_band
354
355ku = maxval(lu_band)        ! largeur de bande sur-diagonale
356kl = maxval(ll_band)        ! largeur de bande sous-diagonale
357
358lbande = kl+ku+1            ! largeur de bande
359ldb    = 2*kl+ku+1          ! pour le dimensionnement de mmat
360
361! eventuellement pour imposer les valeurs de kl et ku
362!kl = klmax
363!ku = kumax
364
365lbande = kl+ku+1            ! largeur de bande
366ldb    = 2*kl+ku+1          ! pour le dimensionnement de mmat
367
368
369! ecriture de la largeur de bande pour en faire eventuellement le trace
370!----------------------------------------------------------------------
371!open(124,file = 'largeur-bande')
372
373!do l = 1,nblig
374!   write(124,*) l,ll_band(l),lu_band(l)
375!end do
376
377!close(124)
378
379! si la matrice mmat  change de taille on la re allocate.
380
381if ((ldb.ne.ldb_old).and.(nblig.ne.nblig_old)) call redim_matrice
382
383
384ldb_old=ldb                 ! pour garder en memoire les anciennes valeurs
385nblig_old=nblig
386
387! ecriture de la matrice proprement dite et du vecteur bdr
388!------------------------------------------------------------
389!$OMP PARALLEL
390!$OMP WORKSHARE
391mmat(:,:) = 0.
392bdr(:,:) = 0.
393!$OMP END WORKSHARE
394
395!$OMP DO PRIVATE(ligne,ilmin,ilmax,jlmin,jlmax,colonne)
396ij_loop : do j = ny1,ny2         ! attention, on rempli 2 lignes par 2 lignes (u,v)
397   do i = nx1,nx2
398
399u_line: if (ok_umat(i,j)) then
400
401   ligne = ligu_L2(i,j)
402
403   if (ligne.lt.0) then          ! teste certains problemes de numeros de ligne
404      write(6,*) 'U_line,i,j',i,j,ligne
405   end if
406
407   bdr(ligne,1) = opposx(i,j)
408
409
410! boucle U : les definitions de ilmin et jlmin sont differentes dans les boucles u et v
411
412      ilmin = max(-1,nx1+1-i)    ! pour avoir i+il entre nx1+1 et nx2
413      ilmax = min(1,nx2-i)
414
415      jlmin = max(-1,ny1-j)      ! pour avoir j+jl entre ny1 et ny2
416      jlmax = min(1,ny2-j)
417
418
419
420iljl_lineU:      do jl = jlmin , jlmax      ! balaye tous les coefficients d'un noeud U de -1 a 1
421         do il = ilmin , ilmax              !                                           de -1 a 1
422
423
424! ecriture des Tu
425               colonne = ligu_L2(i+il,j+jl)
426
427               if (abs(Tu(i,j,il,jl)).gt.eps_col) then
428
429                  if (colonne.lt.0) then    ! trace si un noeud ghost est appele par erreur
430                     write(6,*) 'U_line,i,j,il,jl,Tu',i,j,il,jl,Tu(i,j,il,jl),ghost_x(i+il,j+jl),ghost_y(i+il,j+jl)
431                  end if
432
433
434                  mmat(lbande+ligne-colonne,colonne) = Tu(i,j,il,jl)        ! ecriture des Tu
435
436               endif
437
438! ecriture des Tv
439               colonne = ligv_L2(i+il,j+jl)
440
441               if (abs(Tv(i,j,il,jl)).gt.eps_col) then   
442
443                  if (colonne.lt.0) then  ! trace si un noeud ghost est appele par erreur
444                     write(6,*) 'U_line,i,j,il,jl,Tv',i,j,il,jl,Tv(i,j,il,jl),ghost_x(i+il,j+jl),ghost_y(i+il,j+jl)
445                  end if
446   
447                  mmat(lbande+ligne-colonne,colonne) = Tv(i,j,il,jl)       ! ecriture des Tv
448
449               end if
450
451            end do
452         end do iljl_lineU
453
454      end if u_line
455
456v_line:  if (ok_vmat(i,j)) then
457
458   ligne = ligv_L2(i,j)
459   bdr(ligne,1) = opposy(i,j)
460
461
462! boucle V : les definitions de ilmin et jlmin sont differentes dans les boucles u et v
463
464      ilmin = max(-1,nx1-i)        ! pour avoir i+il entre nx1 et nx2
465      ilmax = min(1,nx2-i)
466
467      jlmin = max(-1,ny1+1-j)      ! pour avoir j+jl entre ny1+1 et ny2
468      jlmax = min(1,ny2-j)
469
470
471iljl_lineV:      do jl = jlmin , jlmax      ! balaye tous les coefficients d'un noeud U de -1 a 1
472         do il = ilmin , ilmax              !                                           de -1 a 1
473
474! ecriture des Su   
475               colonne = ligu_L2(i+il,j+jl)
476
477               if (abs(Su(i,j,il,jl)).gt.eps_col) then
478
479                  if (colonne.lt.0) then  ! trace si un noeud ghost est appele par erreur
480                     write(6,*) 'V_line,i,j,il,jl,Su',i,j,il,jl,Su(i,j,il,jl),ghost_x(i+il,j+jl),ghost_y(i+il,j+jl)
481                  end if
482
483                  mmat(lbande+ligne-colonne,colonne) = Su(i,j,il,jl)      ! ecriture des Su
484             
485               end if
486
487! ecriture des Sv
488               colonne = ligv_L2(i+il,j+jl)
489
490               if (abs(Sv(i,j,il,jl)).gt.eps_col) then
491
492                  if (colonne.lt.0) then  ! trace si un noeud ghost est appele par erreur
493                     write(6,*) 'V_line,i,j,il,jl,Sv',i,j,il,jl,Sv(i,j,il,jl),ghost_x(i+il,j+jl),ghost_y(i+il,j+jl)
494                  end if
495
496
497                  mmat(lbande+ligne-colonne,colonne) = Sv(i,j,il,jl)      ! ecriture des Sv
498
499               end if
500
501            end do
502      end do iljl_lineV
503   end if v_line
504 end do
505end do ij_loop
506!$OMP END DO
507!$OMP END PARALLEL
508
509
510return
511end subroutine prepare_mat
512
513!-----------------------------------------------------------------------------
514subroutine erreur_sgbsv(ifail_sgbsv)
515
516! donne quelques infos en cas d'erreur sgbsv
517
518use module3D_phy,     only: nx,ny
519use remplimat_declar, only: ifail_sgbsv,pos_ligu,pos_ligv,tu,tv,ligu_l2,ligv_l2,il,jl, &
520                            ok_umat,ok_vmat,su,sv
521
522implicit none
523
524integer  :: ifail_sgbsv
525integer  :: li,lj      ! coordonnees du point a probleme
526integer  :: i,j
527
528if (ifail_sgbsv.gt.0) then
529
530   write(*,*) ' ---------------------------------------'
531   write(*,*) '  ERREUR DANS L''ELIMINATION GAUSSIENNE '
532   write(*,*) '  ERREUR Numero',ifail_sgbsv
533   write(*,*) 
534   write(*,*) ' vecteurs colineaires                   '
535   write(*,*) ' diagonale nulle dans decomposition LU  '
536
537! recuperation du noeud de l'erreur
538
539   if (pos_ligu(ifail_sgbsv,1).gt.0) then                              ! pb sur une vitesse U
540      write(*,*) ' vitesse U        i =', pos_ligu(ifail_sgbsv,1), &
541                       '            j =', pos_ligu(ifail_sgbsv,2)
542
543 
544      li=pos_ligu(ifail_sgbsv,1)
545      lj=pos_ligu(ifail_sgbsv,2)
546
547! balaye les noeuds autour  pb avec le format ci-dessous
548997 format(3(i0,2x),'a4',2x,es12.2,2x,L3)
549
550      do j=max(lj-1,1),min(lj+1,ny)
551         do i=max(li-1,1),min(li+1,ny)
552
553            write(6,*)
554            write(6,*)'Tu, Tv, i,j,' ,i,j
555            write(6,*)'--------------------'
556
557               do jl = -1 , 1      ! balaye tous les coefficients de la colonne Uij de -1 a 1
558                  do il = -1 , 1   !                                                de -1 a 1
559 
560                     if ((Tu(i,j,il,jl)).ne.0.) then
561                        write(6,997) il,jl,ligu_L2(i+il,j+jl),'Tu=',Tu(i,j,il,jl), ok_umat(i+il,j+jl)
562                     end if
563                     if ((Tv(i,j,il,jl)).ne.0.) then
564                        write(6,997) il,jl,ligv_L2(i+il,j+jl),'Tv=',Tv(i,j,il,jl), ok_vmat(i+il,j+jl)
565                     end if
566                  end do
567               end do
568
569            end do
570         end do
571
572   else
573      write(*,*) ' vitesse V        i =', pos_ligv(ifail_sgbsv,1), &
574                       '            j =', pos_ligv(ifail_sgbsv,2)
575
576
577      li=pos_ligv(ifail_sgbsv,1)
578      lj=pos_ligv(ifail_sgbsv,2)
579
580      do j=max(lj-1,1),min(lj+1,ny)
581         do i=max(li-1,1),min(li+1,ny)
582           
583            write(6,*)
584            write(6,*)'Su, Sv, i,j,' ,i,j
585
586               do jl = -1 , 1       ! balaye tous les coefficients de la colonne Uij de -1 a 1
587                  do il = -1 , 1    !                                                de -1 a 1
588                     if ((Su(i,j,il,jl)).ne.0.) then
589                       write(6,*) il,jl,ligu_L2(i+il,j+jl),'Su=',Su(i,j,il,jl), ok_umat(i+il,j+jl)
590!                        write(6,997) il,jl,ligu_L2(i+il,j+jl),'Su=',Su(i,j,il,jl), ok_umat(i+il,j+jl)
591
592                     end if
593                     if ((Sv(i,j,il,jl)).ne.0.) then
594                       write(6,*) il,jl,ligv_L2(i+il,j+jl),'Sv=',Sv(i,j,il,jl), ok_vmat(i+il,j+jl)
595!                        write(6,997) il,jl,ligv_L2(i+il,j+jl),'Sv=',Sv(i,j,il,jl), ok_vmat(i+il,j+jl)
596                     end if
597
598                  end do
599               end do
600            end do
601         end do
602
603   endif
604
605   write(*,*) 
606   write(*,*) ' ---------------------------------------'
607 
608
609else if (ifail_sgbsv.lt.0) then
610
611   write(*,*) ' ---------------------------------------'
612   write(*,*) '  ERREUR DANS L''ELIMINATION GAUSSIENNE '
613   write(*,*) '  ERREUR Numero',ifail_sgbsv
614   write(*,*) 
615   write(*,*) ' argument a une valeur illegale         '
616   write(*,*) 
617   write(*,*) ' ---------------------------------------'
618   
619     
620end if
621end subroutine erreur_sgbsv
622!-----------------------------------------------------------------------------
623
624
625end module eq_ellip_sgbsv_mod
626
627!-----------------------------------------------------------------------------
628!
629! rappel arguments sgbsv
630!
631!      SUBROUTINE SGBSV( N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO )
632!
633!  Arguments
634!  =========
635!
636!  N       (input) INTEGER
637!          The number of linear equations, i.e., the order of the
638!          matrix A.  N >= 0.
639!
640!  KL      (input) INTEGER
641!          The number of subdiagonals within the band of A.  KL >= 0.
642!
643!  KU      (input) INTEGER
644!          The number of superdiagonals within the band of A.  KU >= 0.
645!
646!  NRHS    (input) INTEGER
647!          The number of right hand sides, i.e., the number of columns
648!          of the matrix B.  NRHS >= 0.
649!
650!  AB      (input/output) REAL array, dimension (LDAB,N)
651!          On entry, the matrix A in band storage, in rows KL+1 to
652!          2*KL+KU+1; rows 1 to KL of the array need not be set.
653!          The j-th column of A is stored in the j-th column of the
654!          array AB as follows:
655!          AB(KL+KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+KL)
656!          On exit, details of the factorization: U is stored as an
657!          upper triangular band matrix with KL+KU superdiagonals in
658!          rows 1 to KL+KU+1, and the multipliers used during the
659!          factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
660!          See below for further details.
661!
662!  LDAB    (input) INTEGER
663!          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
664!
665!  IPIV    (output) INTEGER array, dimension (N)
666!          The pivot indices that define the permutation matrix P;
667!          row i of the matrix was interchanged with row IPIV(i).
668!
669!  B       (input/output) REAL array, dimension (LDB,NRHS)
670!          On entry, the N-by-NRHS right hand side matrix B.
671!          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
672!
673!  LDB     (input) INTEGER
674!          The leading dimension of the array B.  LDB >= max(1,N).
675!
676!  INFO    (output) INTEGER
677!          = 0:  successful exit
678!          < 0:  if INFO = -i, the i-th argument had an illegal value
679!          > 0:  if INFO = i, U(i,i) is exactly zero.  The factorization
680!                has been completed, but the factor U is exactly
681!                singular, and the solution has not been computed.
682!-----------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.