1 /*
2 $Id$
3 *
4 *
5 *
6 * This is the fortran code wrapper for Choleski factorization
7 * of a positive definite symmetric matrix A in compact column
8 * storage format.
9 *
10 *
11 * subroutine choleski ( n, colQ, mapQ, iwork, work, info)
12 * integer n, mapQ(*), iwork(*), info
13 * DoublePrecision precision colQ(*), work(*)
14 *
15 * Arguments:
16 *
17 * n = (input/integer) dimension of the matrix
18 *
19 * colQ(*) = (input/output/DoublePrecision precison array)
20 * input matrix ( as described in the data structure
21 * of the manual).
22 * The output is the Choleski factor of the matrix
23 * column distributed according to mapQ(*).
24 *
25 * mapQ(*) = (input/integer array) length n, mapQ(i) = process
26 * id of the processor that holds this column.
27 *
28 * iwork(*) = (scratch/integer) length 2n
29 *
30 * work(*) = (scratch/DoublePrecision precision) length 2n + 1
31 * ( actually nprocsQ + n + 1).
32 *
33 * Work must contain at least bufsiz bytes (see cmbbrf.h)
34 *
35 * info = (output/integer) = 0 ; exited normally
36 *
37 */
38
39 #include <stdio.h>
40
41 #include "globalp.c.h"
42
choleski_(n,colQ,mapQ,iwork,work,info)43 void choleski_( n, colQ, mapQ, iwork, work, info)
44 Integer *n, *mapQ, *iwork, *info;
45 DoublePrecision *colQ, *work;
46 {
47 /*
48 this is the Fortran 77 wrapper for the choleski factorization routine
49 the matrix Q and W are assume to be in a column wrapped format
50
51 */
52
53 Integer icont,ibuz;
54 Integer me, nvecsQ;
55 Integer *iscrat, i, j, linfo;
56
57 DoublePrecision **matQ;
58 DoublePrecision *dscrat;
59
60 extern Integer mxmynd_();
61 extern Integer fil_mapvec_();
62 extern void choleski();
63 extern Integer count_list();
64 extern void xerbla_(), pxerbla2_();
65 extern void g_exit_();
66
67 i = 0;
68 me = mxmynd_();
69
70 if( info == NULL ) {
71 i = -6;
72 xerbla_( "Choleski ", &i);
73 }
74
75 if( n == NULL ) {
76 i = -1;
77 xerbla_( "Choleski ", &i);
78 }
79 else
80 if( *n < 1 ) {
81 i = -1;
82 xerbla_( "Choleski ", &i);
83 }
84 else
85 {
86 iscrat = mapQ;
87 for( j = 0; j < *n; j++ ){
88 if( iscrat == NULL ){
89 i = -3;
90 xerbla_( "Choleski \n", &i);
91 }
92 else
93 iscrat++;
94 }
95 }
96
97 /*
98
99 at this point inputs are minimally acceptable for a given processor
100 check to see if n and mapQ are the same set of processors
101
102 */
103
104 iwork[0] = *n;
105 iscrat = iwork + 1;
106 for( i = 0; i < *n; i++ )
107 *(iscrat++) = mapQ[i];
108
109 linfo = 0;
110 i = *n+1;
111 iscrat = iwork + *n + 1;
112
113 i = i * sizeof(Integer);
114 pxerbla2_( &i, (char *) iwork, mapQ, n, iscrat, &linfo );
115
116 g_exit_( &linfo, "Choleski:Mapping inconsistancies.\n", mapQ, n, iwork, work );
117
118 /*
119 if everything seems to agree
120 */
121
122 me = mxmynd_();
123 matQ =(DoublePrecision **) work;
124 nvecsQ = count_list( me, mapQ, n);
125
126 if( nvecsQ < 1) {
127 *info = -50;
128 return;
129 }
130
131 dscrat =(DoublePrecision *)( work + nvecsQ );
132 icont=0;
133 ibuz=0;
134 for( i = 0; i < *n; i++ )
135 {
136 if( mapQ[i] == me){
137 matQ[icont] = &colQ[ibuz];
138 icont++;
139 ibuz+= *n-i;
140 }
141 }
142
143 iscrat = iwork;
144 choleski( n, matQ, mapQ, iscrat, dscrat, info );
145
146 return;
147 }
148