001: /* ///////////////////////////// P /// L /// A /// S /// M /// A /////////////////////////////// */
002: /* ///                    PLASMA computational routines (version 2.1.0)                      ///
003:  * ///                    Author: Jakub Kurzak                                               ///
004:  * ///                    Release Date: November, 15th 2009                                  ///
005:  * ///                    PLASMA is a software package provided by Univ. of Tennessee,       ///
006:  * ///                    Univ. of California Berkeley and Univ. of Colorado Denver          /// */
007: /* ///////////////////////////////////////////////////////////////////////////////////////////// */
008: #include "common.h"
009: 
010: /* /////////////////////////// P /// U /// R /// P /// O /// S /// E /////////////////////////// */
011: // PLASMA_dpotrf - Computes the Cholesky factorization of a symmetric positive definite
012: // (or Hermitian positive definite in the complex case) matrix A.
013: // The factorization has the form
014: //
015: //   A = U**T * U, if uplo = PlasmaUpper, or
016: //   A = L * L**T, if uplo =  PlasmaLower,
017: //
018: // where U is an upper triangular matrix and L is a lower triangular matrix.
019: 
020: /* ///////////////////// A /// R /// G /// U /// M /// E /// N /// T /// S ///////////////////// */
021: // uplo     PLASMA_enum (IN)
022: //          = PlasmaUpper: Upper triangle of A is stored;
023: //          = PlasmaLower: Lower triangle of A is stored.
024: //
025: // N        int (IN)
026: //          The order of the matrix A. N >= 0.
027: //
028: // A        double* (INOUT)
029: //          On entry, the symmetric positive definite (or Hermitian) matrix A.
030: //          If uplo = PlasmaUpper, the leading N-by-N upper triangular part of A
031: //          contains the upper triangular part of the matrix A, and the strictly lower triangular
032: //          part of A is not referenced.
033: //          If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower
034: //          triangular part of the matrix A, and the strictly upper triangular part of A is not
035: //          referenced.
036: //          On exit, if return value = 0, the factor U or L from the Cholesky factorization
037: //          A = U**T*U or A = L*L**T.
038: //
039: // LDA      int (IN)
040: //          The leading dimension of the array A. LDA >= max(1,N).
041: 
042: /* ///////////// R /// E /// T /// U /// R /// N /////// V /// A /// L /// U /// E ///////////// */
043: //          = 0: successful exit
044: //          < 0: if -i, the i-th argument had an illegal value
045: //          > 0: if i, the leading minor of order i of A is not positive definite, so the
046: //               factorization could not be completed, and the solution has not been computed.
047: 
048: /* //////////////////////////////////// C /// O /// D /// E //////////////////////////////////// */
049: int PLASMA_dpotrf(PLASMA_enum uplo, int N, double *A, int LDA)
050: {
051:     int NB, NT;
052:     int status;
053:     double *Abdl;
054:     plasma_context_t *plasma;
055: 
056:     plasma = plasma_context_self();
057:     if (plasma == NULL) {
058:         plasma_fatal_error("PLASMA_dpotrf", "PLASMA not initialized");
059:         return PLASMA_ERR_NOT_INITIALIZED;
060:     }
061:     /* Check input arguments */
062:     if (uplo != PlasmaUpper && uplo != PlasmaLower) {
063:         plasma_error("PLASMA_dpotrf", "illegal value of uplo");
064:         return -1;
065:     }
066:     if (N < 0) {
067:         plasma_error("PLASMA_dpotrf", "illegal value of N");
068:         return -2;
069:     }
070:     if (LDA < max(1, N)) {
071:         plasma_error("PLASMA_dpotrf", "illegal value of LDA");
072:         return -4;
073:     }
074:     /* Quick return */
075:     if (max(N, 0) == 0)
076:         return PLASMA_SUCCESS;
077: 
078:     /* Tune NB depending on M, N & NRHS; Set NBNBSIZE */
079:     status = plasma_tune(PLASMA_FUNC_DPOSV, N, N, 0);
080:     if (status != PLASMA_SUCCESS) {
081:         plasma_error("PLASMA_dpotrf", "plasma_tune() failed");
082:         return status;
083:     }
084: 
085:     /* Set NT */
086:     NB = PLASMA_NB;
087:     NT = (N%NB==0) ? (N/NB) : (N/NB+1);
088: 
089:     /* Allocate memory for matrices in block layout */
090:     Abdl = (double *)plasma_shared_alloc(plasma, NT*NT*PLASMA_NBNBSIZE, PlasmaRealDouble);
091:     if (Abdl == NULL) {
092:         plasma_error("PLASMA_dpotrf", "plasma_shared_alloc() failed");
093:         return PLASMA_ERR_OUT_OF_RESOURCES;
094:     }
095: 
096:     PLASMA_desc descA = plasma_desc_init(
097:         Abdl, PlasmaRealDouble,
098:         PLASMA_NB, PLASMA_NB, PLASMA_NBNBSIZE,
099:         N, N, 0, 0, N, N);
100: 
101:     plasma_parallel_call_3(plasma_lapack_to_tile,
102:         double*, A,
103:         int, LDA,
104:         PLASMA_desc, descA);
105: 
106:     /* Call the native interface */
107:     status = PLASMA_dpotrf_Tile(uplo, &descA);
108: 
109:     if (status == PLASMA_SUCCESS)
110:         plasma_parallel_call_3(plasma_tile_to_lapack,
111:             PLASMA_desc, descA,
112:             double*, A,
113:             int, LDA);
114: 
115:     plasma_shared_free(plasma, Abdl);
116:     return status;
117: }
118: 
119: /* /////////////////////////// P /// U /// R /// P /// O /// S /// E /////////////////////////// */
120: // PLASMA_dpotrf_Tile - Computes the Cholesky factorization of a symmetric positive definite
121: // (or Hermitian positive definite in the complex case) matrix A.
122: // The factorization has the form
123: //
124: //   A = U**T * U, if uplo = PlasmaUpper, or
125: //   A = L * L**T, if uplo =  PlasmaLower,
126: //
127: // where U is an upper triangular matrix and L is a lower triangular matrix.
128: // All matrices are passed through descriptors. All dimensions are taken from the descriptors.
129: 
130: /* ///////////////////// A /// R /// G /// U /// M /// E /// N /// T /// S ///////////////////// */
131: // uplo     PLASMA_enum (IN)
132: //          = PlasmaUpper: Upper triangle of A is stored;
133: //          = PlasmaLower: Lower triangle of A is stored.
134: //
135: // A        double* (INOUT)
136: //          On entry, the symmetric positive definite (or Hermitian) matrix A.
137: //          If uplo = PlasmaUpper, the leading N-by-N upper triangular part of A
138: //          contains the upper triangular part of the matrix A, and the strictly lower triangular
139: //          part of A is not referenced.
140: //          If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower
141: //          triangular part of the matrix A, and the strictly upper triangular part of A is not
142: //          referenced.
143: //          On exit, if return value = 0, the factor U or L from the Cholesky factorization
144: //          A = U**T*U or A = L*L**T.
145: 
146: /* ///////////// R /// E /// T /// U /// R /// N /////// V /// A /// L /// U /// E ///////////// */
147: //          = 0: successful exit
148: //          > 0: if i, the leading minor of order i of A is not positive definite, so the
149: //               factorization could not be completed, and the solution has not been computed.
150: 
151: /* //////////////////////////////////// C /// O /// D /// E //////////////////////////////////// */
152: int PLASMA_dpotrf_Tile(PLASMA_enum uplo, PLASMA_desc *A)
153: {
154:     PLASMA_desc descA = *A;
155:     plasma_context_t *plasma;
156: 
157:     plasma = plasma_context_self();
158:     if (plasma == NULL) {
159:         plasma_fatal_error("PLASMA_dpotrf_Tile", "PLASMA not initialized");
160:         return PLASMA_ERR_NOT_INITIALIZED;
161:     }
162:     /* Check descriptors for correctness */
163:     if (plasma_desc_check(&descA) != PLASMA_SUCCESS) {
164:         plasma_error("PLASMA_dpotrf_Tile", "invalid descriptor");
165:         return PLASMA_ERR_ILLEGAL_VALUE;
166:     }
167:     /* Check input arguments */
168:     if (descA.nb != descA.mb) {
169:         plasma_error("PLASMA_dpotrf_Tile", "only square tiles supported");
170:         return PLASMA_ERR_ILLEGAL_VALUE;
171:     }
172:     if (uplo != PlasmaUpper && uplo != PlasmaLower) {
173:         plasma_error("PLASMA_dpotrf_Tile", "illegal value of uplo");
174:         return -1;
175:     }
176:     /* Quick return */
177: /*
178:     if (max(N, 0) == 0)
179:         return PLASMA_SUCCESS;
180: */
181:     plasma_parallel_call_2(plasma_pdpotrf,
182:         PLASMA_enum, uplo,
183:         PLASMA_desc, descA);
184: 
185:     return PLASMA_INFO;
186: }
187: