SCM

SCM Repository

[matrix] Diff of /pkg/src/cscBlocked.c
ViewVC logotype

Diff of /pkg/src/cscBlocked.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 338, Fri Nov 12 21:17:36 2004 UTC revision 339, Fri Nov 12 21:19:19 2004 UTC
# Line 138  Line 138 
138          error ("Structure of A and A-inverse does not agree");          error ("Structure of A and A-inverse does not agree");
139      }      }
140  }  }
141    
142    /**
143     * Search for the element in a compressed sparse matrix at a given row and column
144     *
145     * @param row row index
146     * @param col column index
147     * @param cp column pointers
148     * @param ci row indices
149     *
150     * @return index of element in ci, if it exists, else -1
151     */
152    static R_INLINE
153    int Ind(int row, int col, const int cp[], const int ci[])
154    {
155        int i, i2 = cp[col + 1];
156        for (i = cp[col]; i < i2; i++)
157            if (ci[i] == row) return i;
158        return -1;
159    }
160    
161    /**
162     * Perform one of the matrix operations
163     *  C := alpha*A*A' + beta*C,
164     * or
165     *  C := alpha*A'*A + beta*C,
166     * where A is a compressed, sparse, blocked matrix and
167     * C is a compressed, sparse, symmetric blocked matrix.
168     *
169     * @param uplo 'U' or 'u' for upper triangular storage, else lower.
170     * @param trans 'T' or 't' for transpose.
171     * @param alpha scalar multiplier of outer product
172     * @param A compressed sparse blocked matrix
173     * @param beta scalar multiplier of c
174     * @param C compressed sparse blocked symmetric matrix to be updated
175     */
176    SEXP cscBlocked_syrk(SEXP uplo, SEXP trans, SEXP alpha, SEXP A, SEXP beta, SEXP C)
177    {
178        SEXP axp = GET_SLOT(A, Matrix_xSym),
179            app = GET_SLOT(A, Matrix_pSym),
180            cxp = GET_SLOT(C, Matrix_xSym),
181            cpp = GET_SLOT(C, Matrix_pSym);
182        int *adims = INTEGER(getAttrib(axp, R_DimSymbol)),
183            *ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
184            *ap = INTEGER(app),
185            *cdims = INTEGER(getAttrib(cxp, R_DimSymbol)),
186            *ci = INTEGER(GET_SLOT(A, Matrix_iSym)),
187            *cp = INTEGER(cpp),
188            j, k,
189            nca = length(app) - 1,
190            ncc = length(cpp) - 1,
191            npairs,
192            inplace;
193        char *ul, *tr;
194        double *ax = REAL(axp),
195            *cx = REAL(cxp),
196            bta = asReal(beta);
197    
198        if (length(uplo) != 1 || length(trans) != 1 || length(alpha) != 1 ||
199            length(beta) != 1 || !isReal(alpha) || !isReal(beta) ||
200            !isString(uplo) || !isString(trans))
201            error("uplo and trans must be character scalars, alpha and beta real scalars");
202        if (cdims[0] != cdims[1]) error("blocks in C must be square");
203        ul = CHAR(STRING_ELT(uplo, 0));
204        tr = CHAR(STRING_ELT(trans, 0));
205        if (toupper(tr[0]) == 'T')
206            error("Code for trans == 'T' not yet written");
207    /* FIXME: Write the other version */
208        if (adims[0] != cdims[0])
209            error("Dimension inconsistency in blocks: dim(A)=[%d,,], dim(C)=[%d,,]",
210                  adims[0], cdims[0]);
211                                    /* check the row indices */
212        for (k = 0; k < adims[2]; k++) {
213            int aik = ai[k];
214            if (aik < 0 || aik >= ncc)
215                error("Row index %d = %d is out of range [0, %d]",
216                      k, ai[k], ncc - 1);
217        }
218                                    /* check if C can be updated in place */
219        inplace = 1;
220        npairs = 0;
221        for (j = 0; j < nca; j++) {
222            int k, kk, k2 = ap[j+1];
223            int nnz = k2 - ap[j];
224            npairs += (nnz * (nnz - 1)) / 2;
225            for (k = ap[j]; k < k2; k++) {
226    /* FIXME: This check assumes uplo == 'L' */
227                for (kk = k; kk < k2; kk++) {
228                    if (Ind(ai[k], ai[kk], cp, ci) < 0) {
229                        inplace = 0;
230                        break;
231                    }
232                }
233                if (!inplace) break;
234            }
235        }
236                                    /* multiply C by beta */
237        for (j = 0; j < cdims[0]*cdims[1]*cdims[2]; j++) cx[j] *= bta;
238        if (inplace) {
239            int scalar = (adims[0] == 1 && adims[1] == 1);
240    
241            for (j = 0; j < nca; j++) {
242                int k, kk, k2 = ap[j+1];
243                for (k = ap[j]; k < k2; k++) {
244                    int ii = ai[k];
245                    for (kk = k; kk < k2; kk++) {
246                        int jj = ai[kk], K = Ind(ii, jj, cp, ci);
247                        if (scalar) cx[K] += ax[k] * ax[kk];
248                        else {
249                        }
250                    }
251                }
252            }
253        }
254    }
255    

Legend:
Removed from v.338  
changed lines
  Added in v.339

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge