source: trunk/SOURCES/New-remplimat/eq_ellipt_sgbsv_mod-0.2.f90 @ 4

Last change on this file since 4 was 4, checked in by dumas, 10 years ago

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